From 2ce15a9b1ddbe7832995027c78c5c96837c7420f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 30 Nov 2022 16:00:00 -0500 Subject: [PATCH] Removed most of the rest of the flat binary file cruft --- .../Basic_Simulation/initial_conditions.ipynb | 8 +- .../Basic_Simulation/initial_conditions.py | 6 +- src/CMakeLists.txt | 2 - src/discard/discard.f90 | 4 +- src/encounter/encounter_io.f90 | 102 --- src/io/io.f90 | 701 ++---------------- src/main/swiftest_driver.f90 | 6 +- src/modules/encounter_classes.f90 | 10 - src/modules/rmvs_classes.f90 | 8 - src/modules/swiftest_classes.f90 | 105 +-- src/modules/swiftest_globals.f90 | 10 - src/netcdf/netcdf.f90 | 34 +- src/rmvs/rmvs_io.f90 | 56 -- src/rmvs/rmvs_step.f90 | 12 +- src/setup/setup.f90 | 4 +- src/symba/symba_discard.f90 | 2 +- src/symba/symba_io.f90 | 69 +- src/whm/whm_setup.f90 | 2 +- 18 files changed, 99 insertions(+), 1042 deletions(-) delete mode 100644 src/encounter/encounter_io.f90 delete mode 100644 src/rmvs/rmvs_io.f90 diff --git a/examples/Basic_Simulation/initial_conditions.ipynb b/examples/Basic_Simulation/initial_conditions.ipynb index 2bcd9cfe7..cdcaae0ed 100644 --- a/examples/Basic_Simulation/initial_conditions.ipynb +++ b/examples/Basic_Simulation/initial_conditions.ipynb @@ -21,7 +21,7 @@ "outputs": [], "source": [ "# Initialize the simulation object as a variable\n", - "sim = swiftest.Simulation(tstart=0.0, tstop=1.0e6, dt=0.005, tstep_out=1.0e3, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)" + "sim = swiftest.Simulation(tstart=0.0, tstop=1.0e6, dt=0.01, tstep_out=1.0e3, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)" ] }, { @@ -79,7 +79,7 @@ "metadata": {}, "outputs": [], "source": [ - "sim.add_body(name_pl, a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl, GMpl=GM_pl, Rpl=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl)" + "sim.add_body(name_pl, a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl, Gmass=GM_pl, radius=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl)" ] }, { @@ -144,9 +144,9 @@ ], "metadata": { "kernelspec": { - "display_name": "swiftest", + "display_name": "Python (My debug_env Kernel)", "language": "python", - "name": "swiftest" + "name": "debug_env" }, "language_info": { "codemirror_mode": { diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 4accb67e7..9bb279b51 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -23,7 +23,7 @@ from numpy.random import default_rng # Initialize the simulation object as a variable -sim = swiftest.Simulation(tstart=0.0, tstop=1.0e3, dt=0.005, tstep_out=1.0e0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) +sim = swiftest.Simulation(tstart=0.0, tstop=1.0e3, dt=0.010, tstep_out=1.0e0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) # Add the modern planets and the Sun using the JPL Horizons Database sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"]) @@ -49,7 +49,7 @@ roty_pl = [0.0, 0.0, 0.0, 0.0, 0.0] rotz_pl = [0.0, 0.0, 0.0, 0.0, 0.0] -sim.add_body(name_pl, a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl, GMpl=GM_pl, Rpl=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl) +sim.add_body(name=name_pl, v1=a_pl, v2=e_pl, v3=inc_pl, v4=capom_pl, v5=omega_pl, v6=capm_pl, Gmass=GM_pl, radius=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl) # Add 10 user-defined test particles ntp = 10 @@ -62,7 +62,7 @@ omega_tp = default_rng().uniform(0.0, 360.0, ntp) capm_tp = default_rng().uniform(0.0, 360.0, ntp) -sim.add_body(name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp) +sim.add_body(name=name_tp, v1=a_tp, v2=e_tp, v3=inc_tp, v4=capom_tp, v5=omega_tp, v6=capm_tp) # Display the run configuration parameters sim.get_parameter() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fd0bf10c3..63c89f2b3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -28,7 +28,6 @@ SET(FAST_MATH_FILES ${SRC}/discard/discard.f90 ${SRC}/drift/drift.f90 ${SRC}/encounter/encounter_check.f90 - ${SRC}/encounter/encounter_io.f90 ${SRC}/encounter/encounter_setup.f90 ${SRC}/encounter/encounter_util.f90 ${SRC}/fraggle/fraggle_generate.f90 @@ -53,7 +52,6 @@ SET(FAST_MATH_FILES ${SRC}/orbel/orbel.f90 ${SRC}/rmvs/rmvs_discard.f90 ${SRC}/rmvs/rmvs_encounter_check.f90 - ${SRC}/rmvs/rmvs_io.f90 ${SRC}/rmvs/rmvs_setup.f90 ${SRC}/rmvs/rmvs_step.f90 ${SRC}/rmvs/rmvs_util.f90 diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index ffaba5dd6..6c799844f 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -38,8 +38,8 @@ module subroutine discard_system(self, param) call tp%discard(system, param) ltp_discards = (tp_discards%nbody > 0) end if - - if (lpl_discards .or. ltp_discards) call system%write_discard(param) + if (ltp_discards) call tp_discards%write_particle_info(param%nciu, param) + if (lpl_discards) call pl_discards%write_particle_info(param%nciu, param) if (lpl_discards .and. param%lenergy) call self%conservation_report(param, lterminal=.false.) if (lpl_check) call pl_discards%setup(0,param) if (ltp_check) call tp_discards%setup(0,param) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 deleted file mode 100644 index 199169b35..000000000 --- a/src/encounter/encounter_io.f90 +++ /dev/null @@ -1,102 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (encounter_classes) s_encounter_io - use swiftest -contains - - - module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) - !! author: David A. Minton - !! - !! Write a single frame of close encounter data to output binary files - !! - !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 - !! Adapted from Hal Levison's Swift routine io_write_encounter.f - implicit none - ! Arguments - integer(I4B), intent(in) :: iu !! Open file unit number - real(DP), intent(in) :: t !! Time of encounter - integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies - real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies - real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies - real(DP), dimension(:), intent(in) :: xh1, xh2 !! Heliocentric position vectors of the two encountering bodies - real(DP), dimension(:), intent(in) :: vh1, vh2 !! Heliocentric velocity vectors of the two encountering bodies - ! Internals - character(len=STRMAX) :: errmsg - - write(iu, err=667, iomsg=errmsg) t - write(iu, err=667, iomsg=errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 - write(iu, err=667, iomsg=errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 - - return - 667 continue - write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine - - module subroutine encounter_io_write_list(self, pl, encbody, param) - !! author: David A. Minton - !! - !! Write the encounters to the output encounter binary files - implicit none - ! Arguments - class(encounter_list), intent(in) :: self !! Swiftest encounter list object - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - logical , save :: lfirst = .true. - integer(I4B) :: k, ierr - character(len=STRMAX) :: errmsg - - if (param%enc_out == "" .or. self%nenc == 0) return - - open(unit=LUN, file=param%enc_out, status='OLD', position='APPEND', form='UNFORMATTED', iostat=ierr, iomsg=errmsg) - if (ierr /= 0) then - if (lfirst) then - open(unit=LUN, file=param%enc_out, status='NEW', form='UNFORMATTED', err=667, iomsg=errmsg) - else - goto 667 - end if - end if - lfirst = .false. - - associate(ind1 => self%index1, ind2 => self%index2) - select type(encbody) - class is (swiftest_pl) - do k = 1, self%nenc - call encounter_io_write_frame(LUN, self%t(k), & - pl%id(ind1(k)), encbody%id(ind2(k)), & - pl%Gmass(ind1(k)), encbody%Gmass(ind2(k)), & - pl%radius(ind1(k)), encbody%radius(ind2(k)), & - self%x1(:,k), self%x2(:,k), & - self%v1(:,k), self%v2(:,k)) - end do - class is (swiftest_tp) - do k = 1, self%nenc - call encounter_io_write_frame(LUN, self%t(k), & - pl%id(ind1(k)), encbody%id(ind2(k)), & - pl%Gmass(ind1(k)), 0.0_DP, & - pl%radius(ind1(k)), 0.0_DP, & - self%x1(:,k), self%x2(:,k), & - self%v1(:,k), self%v2(:,k)) - end do - end select - end associate - - close(unit = LUN, err = 667, iomsg = errmsg) - - return - 667 continue - write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine encounter_io_write_list - -end submodule s_encounter_io \ No newline at end of file diff --git a/src/io/io.f90 b/src/io/io.f90 index 2d7f59de2..de4c547e9 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -179,10 +179,7 @@ module subroutine io_conservation_report(self, param, lterminal) write(*,*) "Severe error! Mass not conserved! Halting!" ! Save the frame of data to the bin file in the slot just after the present one for diagnostics param%ioutput = param%ioutput + 1_I8B - call pl%xv2el(cb) - call self%write_hdr(param%nciu, param) - call cb%write_frame(param%nciu, param) - call pl%write_frame(param%nciu, param) + call self%write_frame(param%nciu, param) call param%nciu%close() call util_exit(FAILURE) end if @@ -228,48 +225,6 @@ module subroutine io_dump_param(self, param_file_name) end subroutine io_dump_param - module subroutine io_dump_base(self, param) - !! author: David A. Minton - !! - !! Dump massive body data to files - !! - !! Adapted from David E. Kaufmann's Swifter routine: io_dump_pl.f90 and io_dump_tp.f90 - !! Adapted from Hal Levison's Swift routine io_dump_pl.f and io_dump_tp.f - implicit none - ! Arguments - class(swiftest_base), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: iu = LUN - character(len=:), allocatable :: dump_file_name - character(STRMAX) :: errmsg - - select type(self) - class is(swiftest_cb) - dump_file_name = trim(adjustl(param%incbfile)) - class is (swiftest_pl) - dump_file_name = trim(adjustl(param%inplfile)) - class is (swiftest_tp) - dump_file_name = trim(adjustl(param%intpfile)) - end select - open(unit = iu, file = dump_file_name, form = "UNFORMATTED", status = 'replace', err = 667, iomsg = errmsg) - select type(self) - class is (swiftest_body) - write(iu, err = 667, iomsg = errmsg) self%nbody - call io_write_frame_body(self,iu, param) - class is (swiftest_cb) - call io_write_frame_cb(self,iu, param) - end select - close(iu, err = 667, iomsg = errmsg) - - return - - 667 continue - write(*,*) "Error dumping body data to file " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_dump_base - - module subroutine io_dump_system(self, param) !! author: David A. Minton !! @@ -289,9 +244,9 @@ module subroutine io_dump_system(self, param) allocate(dump_param, source=param) param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) - dump_param%in_form = XV + dump_param%in_form = "XV" dump_param%out_stat = 'APPEND' - dump_param%in_type = NETCDF_DOUBLE_TYPE + dump_param%in_type = "NETCDF_DOUBLE" dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody dump_param%nciu%time_chunk = 1 @@ -299,14 +254,11 @@ module subroutine io_dump_system(self, param) call dump_param%dump(param_file_name) - dump_param%out_form = XV + dump_param%out_form = "XV" dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx))) dump_param%ioutput = 0 call dump_param%nciu%initialize(dump_param) - call self%write_hdr(dump_param%nciu, dump_param) - call self%cb%write_frame(dump_param%nciu, dump_param) - call self%pl%write_frame(dump_param%nciu, dump_param) - call self%tp%write_frame(dump_param%nciu, dump_param) + call self%write_frame(dump_param%nciu, dump_param) call dump_param%nciu%close() ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) call param%nciu%flush(param) @@ -392,45 +344,6 @@ module subroutine io_get_args(integrator, param_file_name, display_style) end subroutine io_get_args - module function io_get_old_t_final_system(self, param) result(old_t_final) - !! author: David A. Minton - !! - !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. - !! - implicit none - ! Arguments - class(swiftest_nbody_system), intent(in) :: self - class(swiftest_parameters), intent(in) :: param - ! Result - real(DP) :: old_t_final - ! Internals - class(swiftest_nbody_system), allocatable :: tmpsys - class(swiftest_parameters), allocatable :: tmpparam - integer(I4B) :: ierr, iu = LUN - character(len=STRMAX) :: errmsg - - old_t_final = 0.0_DP - allocate(tmpsys, source=self) - allocate(tmpparam, source=param) - - ierr = 0 - open(unit = iu, file = param%outfile, status = 'OLD', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - do - ierr = tmpsys%read_frame(iu, tmpparam) - if (ierr /= 0) exit - end do - if (is_iostat_end(ierr)) then - old_t_final = tmpparam%t - close(iu) - return - end if - - 667 continue - write(*,*) "Error reading binary output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end function io_get_old_t_final_system - - module function io_get_token(buffer, ifirst, ilast, ierr) result(token) !! author: David A. Minton !! @@ -754,8 +667,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - if ((param%in_type /= REAL8_TYPE) .and. (param%in_type /= "ASCII") & - .and. (param%in_type /= NETCDF_FLOAT_TYPE) .and. (param%in_type /= NETCDF_DOUBLE_TYPE)) then + if ((param%in_type /= "ASCII") .and. (param%in_type /= "NETCDF_FLOAT") .and. (param%in_type /= "NETCDF_DOUBLE")) then write(iomsg,*) 'Invalid input file type:',trim(adjustl(param%in_type)) iostat = -1 return @@ -772,7 +684,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end if param%lrestart = (param%out_stat == "APPEND") if (param%outfile /= "") then - if ((param%out_type /= NETCDF_FLOAT_TYPE) .and. (param%out_type /= NETCDF_DOUBLE_TYPE)) then + if ((param%out_type /= "NETCDF_FLOAT") .and. (param%out_type /= "NETCDF_DOUBLE")) then write(iomsg,*) 'Invalid out_type: ',trim(adjustl(param%out_type)) iostat = -1 return @@ -949,11 +861,11 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) call io_param_writer_one("TSTOP", param%tstop, unit) call io_param_writer_one("DT", param%dt, unit) call io_param_writer_one("IN_TYPE", param%in_type, unit) - if ((param%in_type == REAL4_TYPE) .or. (param%in_type == REAL8_TYPE)) then + if (param%in_type == "ASCII") then call io_param_writer_one("CB_IN", param%incbfile, unit) call io_param_writer_one("PL_IN", param%inplfile, unit) call io_param_writer_one("TP_IN", param%intpfile, unit) - else if ((param%in_type == NETCDF_FLOAT_TYPE) .or. (param%in_type == NETCDF_DOUBLE_TYPE)) then + else call io_param_writer_one("NC_IN", param%in_netcdf, unit) end if @@ -1204,7 +1116,7 @@ module subroutine io_read_in_base(self,param) class(swiftest_base), intent(inout) :: self !! Swiftest base object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - if ((param%in_type == NETCDF_FLOAT_TYPE) .or. (param%in_type == NETCDF_DOUBLE_TYPE)) return ! This method is not used in NetCDF mode, as reading is done for the whole system, not on individual particle types + if (param%in_type /= "ASCII") return ! This method is not used in NetCDF mode, as reading is done for the whole system, not on individual particle types select type(self) class is (swiftest_body) @@ -1231,13 +1143,13 @@ subroutine io_read_in_body(self, param) ! Internals integer(I4B) :: iu = LUN integer(I4B) :: i, nbody - logical :: is_ascii character(len=:), allocatable :: infile character(STRMAX) :: errmsg ! Internals integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Select the appropriate polymorphic class (test particle or massive body) + if (param%in_type /= "ASCII") return ! Not for NetCDF select type(self) class is (swiftest_pl) @@ -1246,18 +1158,8 @@ subroutine io_read_in_body(self, param) infile = param%intpfile end select - is_ascii = (param%in_type == 'ASCII') - select case(param%in_type) - case(ASCII_TYPE) - open(unit = iu, file = infile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg) - read(iu, *, err = 667, iomsg = errmsg) nbody - case (REAL4_TYPE, REAL8_TYPE) - open(unit=iu, file=infile, status='old', form='UNFORMATTED', err = 667, iomsg = errmsg) - read(iu, err = 667, iomsg = errmsg) nbody - case default - write(errmsg,*) trim(adjustl(param%in_type)) // ' is an unrecognized file type' - goto 667 - end select + open(unit = iu, file = infile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg) + read(iu, *, err = 667, iomsg = errmsg) nbody call self%setup(nbody, param) ierr = 0 @@ -1271,7 +1173,6 @@ subroutine io_read_in_body(self, param) end if close(iu, err = 667, iomsg = errmsg) - if (ierr == 0) return 667 continue @@ -1297,26 +1198,23 @@ subroutine io_read_in_cb(self, param) integer(I4B) :: ierr character(len=NAMELEN) :: name - if (param%in_type == 'ASCII') then - 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) name - call self%info%set_value(name=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 - read(iu, *, err = 667, iomsg = errmsg) self%j2rp2 - read(iu, *, err = 667, iomsg = errmsg) self%j4rp4 - if (param%lrotation) then - read(iu, *, err = 667, iomsg = errmsg) self%Ip(1), self%Ip(2), self%Ip(3) - read(iu, *, err = 667, iomsg = errmsg) self%rot(1), self%rot(2), self%rot(3) - end if - ierr = 0 - else - open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - ierr = self%read_frame(iu, param) + if (param%in_type /= "ASCII") return ! Not for NetCDF + + 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) name + call self%info%set_value(name=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 + read(iu, *, err = 667, iomsg = errmsg) self%j2rp2 + read(iu, *, err = 667, iomsg = errmsg) self%j4rp4 + if (param%lrotation) then + read(iu, *, err = 667, iomsg = errmsg) self%Ip(1), self%Ip(2), self%Ip(3) + read(iu, *, err = 667, iomsg = errmsg) self%rot(1), self%rot(2), self%rot(3) end if + ierr = 0 close(iu, err = 667, iomsg = errmsg) if (ierr == 0) then @@ -1354,18 +1252,7 @@ module subroutine io_read_in_system(self, param) integer(I4B) :: ierr class(swiftest_parameters), allocatable :: tmp_param - if ((param%in_type == NETCDF_DOUBLE_TYPE) .or. (param%in_type == NETCDF_FLOAT_TYPE)) then - allocate(tmp_param, source=param) - tmp_param%outfile = param%in_netcdf - tmp_param%out_form = param%in_form - if (.not. param%lrestart) then - ! Turn off energy computation so we don't have to feed it into the initial conditions - tmp_param%lenergy = .false. - end if - ierr = self%read_frame(tmp_param%nciu, tmp_param) - deallocate(tmp_param) - if (ierr /=0) call util_exit(FAILURE) - else + if (param%in_type == "ASCII") then call self%cb%read_in(param) call self%pl%read_in(param) call self%tp%read_in(param) @@ -1378,6 +1265,17 @@ module subroutine io_read_in_system(self, param) self%Lescape(:) = param%Lescape(:) self%Ecollisions = param%Ecollisions self%Euntracked = param%Euntracked + else + allocate(tmp_param, source=param) + tmp_param%outfile = param%in_netcdf + tmp_param%out_form = param%in_form + if (.not. param%lrestart) then + ! Turn off energy computation so we don't have to feed it into the initial conditions + tmp_param%lenergy = .false. + end if + ierr = self%read_frame(tmp_param%nciu, tmp_param) + deallocate(tmp_param) + if (ierr /=0) call util_exit(FAILURE) end if param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP)) @@ -1390,56 +1288,6 @@ module subroutine io_read_in_system(self, param) end subroutine io_read_in_system - function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & - xh1, xh2, vh1, vh2, enc_out, out_type) result(ierr) - !! author: David A. Minton - !! - !! Read close encounter data from input binary files - !! Other than time t, there is no direct file input from this function - !! Function returns read error status (0 = OK, nonzero = ERROR) - !! Adapted from David E. Kaufmann's Swifter routine: io_read_encounter.f90 - implicit none - ! Arguments - integer(I4B), intent(out) :: id1, id2 - real(DP), intent(out) :: t, Gmass1, Gmass2, radius1, radius2 - real(DP), dimension(:), intent(out) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: enc_out, out_type - ! Result - integer(I4B) :: ierr - ! Internals - logical , save :: lfirst = .true. - integer(I4B), save :: iu = lun - - if (lfirst) then - open(unit = iu, file = enc_out, status = 'OLD', form = 'UNFORMATTED', iostat = ierr) - if (ierr /= 0) then - write(*, *) "Swiftest Error:" - write(*, *) " unable to open binary encounter file" - call util_exit(FAILURE) - end if - lfirst = .false. - end if - read(iu, iostat = ierr) t - if (ierr /= 0) then - close(unit = iu, iostat = ierr) - return - end if - - read(iu, iostat = ierr) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), vh1(3), Gmass1, radius1 - if (ierr /= 0) then - close(unit = iu, iostat = ierr) - return - end if - read(iu, iostat = ierr) id2, xh2(2), xh2(2), xh2(3), vh2(2), vh2(2), vh2(3), Gmass2, radius2 - if (ierr /= 0) then - close(unit = iu, iostat = ierr) - return - end if - - return - end function io_read_encounter - - module function io_read_frame_body(self, iu, param) result(ierr) !! author: David A. Minton !! @@ -1462,14 +1310,14 @@ module function io_read_frame_body(self, iu, param) result(ierr) if (self%nbody == 0) return - if ((param%in_form /= EL) .and. (param%in_form /= XV)) then + if ((param%in_form /= "EL") .and. (param%in_form /= "XV")) then write(errmsg, *) trim(adjustl(param%in_form)) // " is not a recognized format code for input files." goto 667 end if associate(n => self%nbody) - if (param%in_form == EL) then + if (param%in_form == "EL") then if (.not.allocated(self%a)) allocate(self%a(n)) if (.not.allocated(self%e)) allocate(self%e(n)) if (.not.allocated(self%inc)) allocate(self%inc(n)) @@ -1479,53 +1327,7 @@ module function io_read_frame_body(self, iu, param) result(ierr) end if select case(param%in_type) - case (REAL4_TYPE, REAL8_TYPE) - read(iu, err = 667, iomsg = errmsg) self%id(:) - read(iu, err = 667, iomsg = errmsg) name(:) - do i = 1, n - call self%info(i)%set_value(name=name(i)) - end do - - select case (param%in_form) - case (XV) - read(iu, err = 667, iomsg = errmsg) self%xh(1, :) - read(iu, err = 667, iomsg = errmsg) self%xh(2, :) - read(iu, err = 667, iomsg = errmsg) self%xh(3, :) - read(iu, err = 667, iomsg = errmsg) self%vh(1, :) - read(iu, err = 667, iomsg = errmsg) self%vh(2, :) - read(iu, err = 667, iomsg = errmsg) self%vh(3, :) - case (EL) - read(iu, err = 667, iomsg = errmsg) self%a(:) - read(iu, err = 667, iomsg = errmsg) self%e(:) - read(iu, err = 667, iomsg = errmsg) self%inc(:) - read(iu, err = 667, iomsg = errmsg) self%capom(:) - read(iu, err = 667, iomsg = errmsg) self%omega(:) - read(iu, err = 667, iomsg = errmsg) self%capm(:) - end select - - select type(pl => self) - class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body - read(iu, err = 667, iomsg = errmsg) pl%Gmass(:) - pl%mass(:) = pl%Gmass(:) / param%GU - if (param%lrhill_present) read(iu, err = 667, iomsg = errmsg) pl%rhill(:) - if (param%lclose) read(iu, err = 667, iomsg = errmsg) pl%radius(:) - if (param%lrotation) then - read(iu, err = 667, iomsg = errmsg) pl%Ip(1, :) - read(iu, err = 667, iomsg = errmsg) pl%Ip(2, :) - read(iu, err = 667, iomsg = errmsg) pl%Ip(3, :) - read(iu, err = 667, iomsg = errmsg) pl%rot(1, :) - read(iu, err = 667, iomsg = errmsg) pl%rot(2, :) - read(iu, err = 667, iomsg = errmsg) pl%rot(3, :) - end if - ! if (param%ltides) then - ! read(iu, err = 667, iomsg = errmsg) pl%k2(:) - ! read(iu, err = 667, iomsg = errmsg) pl%Q(:) - ! end if - end select - - param%maxid = max(param%maxid, maxval(self%id(1:n))) - - case (ASCII_TYPE) + case ("ASCII") do i = 1, n select type(self) class is (swiftest_pl) @@ -1543,10 +1345,10 @@ module function io_read_frame_body(self, iu, param) result(ierr) call self%info(i)%set_value(name=name(i)) select case(param%in_form) - case (XV) + case ("XV") read(iu, *, err = 667, iomsg = errmsg) self%xh(1, i), self%xh(2, i), self%xh(3, i) read(iu, *, err = 667, iomsg = errmsg) self%vh(1, i), self%vh(2, i), self%vh(3, i) - case (EL) + case ("EL") read(iu, *, err = 667, iomsg = errmsg) self%a(i), self%e(i), self%inc(i) read(iu, *, err = 667, iomsg = errmsg) self%capom(i), self%omega(i), self%capm(i) end select @@ -1567,7 +1369,7 @@ module function io_read_frame_body(self, iu, param) result(ierr) end do end select - if (param%in_form == EL) then + if (param%in_form == "EL") then self%inc(1:n) = self%inc(1:n) * DEG2RAD self%capom(1:n) = self%capom(1:n) * DEG2RAD self%omega(1:n) = self%omega(1:n) * DEG2RAD @@ -1591,146 +1393,6 @@ module function io_read_frame_body(self, iu, param) result(ierr) end function io_read_frame_body - module function io_read_frame_cb(self, iu, param) result(ierr) - !! author: David A. Minton - !! - !! Reads a frame of output of central body data to the binary output file - !! - !! Adapted from David E. Kaufmann's Swifter routine io_read_frame.f90 - !! Adapted from Hal Levison's Swift routine io_read_frame.F - implicit none - ! Arguments - class(swiftest_cb), intent(inout) :: self !! Swiftest central body 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 - ! Result - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful - ! Internals - character(len=STRMAX) :: errmsg - character(len=NAMELEN) :: name - - read(iu, err = 667, iomsg = errmsg) self%id - read(iu, err = 667, iomsg = errmsg) name - call self%info%set_value(name=name) - read(iu, err = 667, iomsg = errmsg) self%Gmass - self%mass = self%Gmass / param%GU - read(iu, err = 667, iomsg = errmsg) self%radius - read(iu, err = 667, iomsg = errmsg) self%j2rp2 - read(iu, err = 667, iomsg = errmsg) self%j4rp4 - if (param%lrotation) then - read(iu, err = 667, iomsg = errmsg) self%Ip(1) - read(iu, err = 667, iomsg = errmsg) self%Ip(2) - read(iu, err = 667, iomsg = errmsg) self%Ip(3) - read(iu, err = 667, iomsg = errmsg) self%rot(1) - read(iu, err = 667, iomsg = errmsg) self%rot(2) - read(iu, err = 667, iomsg = errmsg) self%rot(3) - end if - ! if (param%ltides) then - ! read(iu, err = 667, iomsg = errmsg) self%k2 - ! read(iu, err = 667, iomsg = errmsg) self%Q - ! end if - - ierr = 0 - return - - 667 continue - write(*,*) "Error reading central body frame: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end function io_read_frame_cb - - - module function io_read_frame_system(self, iu, param) result(ierr) - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Read a frame (header plus records for each massive body and active test particle) from a output binary file - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system 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 - ! Result - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful - ! Internals - character(len=STRMAX) :: errmsg - integer(I4B) :: npl, ntp - - ierr = io_read_hdr(iu, param%t, npl, ntp, param%out_form, param%out_type) - if (is_iostat_end(ierr)) return ! Reached the end of the frames - call self%pl%setup(npl, param) - call self%tp%setup(ntp, param) - - if (ierr /= 0) then - write(errmsg, *) "Cannot read frame header." - goto 667 - end if - ierr = self%cb%read_frame(iu, param) - if (ierr /= 0) then - write(errmsg, *) "Cannot read central body frame." - goto 667 - end if - ierr = self%pl%read_frame(iu, param) - if (ierr /= 0) then - write(errmsg, *) "Cannot read massive body frame." - goto 667 - end if - ierr = self%tp%read_frame(iu, param) - if (ierr /= 0) then - write(errmsg, *) "Cannot read test particle frame." - goto 667 - end if - - return - - 667 continue - write(*,*) "Error reading system frame: " // trim(adjustl(errmsg)) - end function io_read_frame_system - - - function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) - !! author: David A. Minton - !! - !! Read frame header from input binary files - !! Function returns read error status (0 = OK, nonzero = ERROR) - !! Adapted from David E. Kaufmann's Swifter routine: io_read_hdr.f90 - !! Adapted from Hal Levison's Swift routine io_read_hdr.f - implicit none - ! Arguments - integer(I4B), intent(in) :: iu - integer(I4B), intent(out) :: npl, ntp - character(*), intent(out) :: out_form - real(DP), intent(out) :: t - character(*), intent(in) :: out_type - ! Result - integer(I4B) :: ierr - ! Internals - real(SP) :: ttmp - character(len=STRMAX) :: errmsg - - select case (out_type) - case (REAL4_TYPE) - read(iu, iostat = ierr, err = 667, iomsg = errmsg, end = 333) ttmp - t = ttmp - case (REAL8_TYPE) - read(iu, iostat = ierr, err = 667, iomsg = errmsg, end = 333) t - case default - write(errmsg,*) trim(adjustl(out_type)) // ' is an unrecognized file type' - ierr = -1 - end select - read(iu, iostat = ierr, err = 667, iomsg = errmsg) npl - read(iu, iostat = ierr, err = 667, iomsg = errmsg) ntp - read(iu, iostat = ierr, err = 667, iomsg = errmsg) out_form - - return - - 667 continue - write(*,*) "Error reading header: " // trim(adjustl(errmsg)) - 333 continue - return - - return - end function io_read_hdr - - module subroutine io_read_in_param(self, param_file_name) !! author: David A. Minton !! @@ -1763,33 +1425,6 @@ 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 - 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%collision_id - read(iu, err = 667, iomsg = errmsg) self%origin_xh(:) - read(iu, err = 667, iomsg = errmsg) self%origin_vh(:) - - 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_set_display_param(self, display_style) !! author: David A. Minton !! @@ -1850,197 +1485,6 @@ module subroutine io_toupper(string) end subroutine io_toupper - module subroutine io_write_discard(self, param) - !! author: David A. Minton - !! - !! Write out information about discarded test particle - !! - !! Adapted from David E. Kaufmann's Swifter routine io_discard_write.f90 - !! Adapted from Hal Levison's Swift routine io_discard_write.f - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i - logical, save :: lfirst = .true. - real(DP), dimension(:,:), allocatable :: vh - character(*), parameter :: HDRFMT = '(E23.16, 1X, I8, 1X, L1)' - character(*), parameter :: NAMEFMT = '(A, 2(1X, I8))' - character(*), parameter :: VECFMT = '(3(E23.16, 1X))' - character(*), parameter :: NPLFMT = '(I8)' - character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' - class(swiftest_body), allocatable :: pltemp - character(len=STRMAX) :: errmsg, out_stat - - associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody) - - ! Record the discarded body metadata information to file - if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - call tp_discards%write_particle_info(param%nciu, param) - end if - - if (param%discard_out == "") return - - if (nsp == 0) return - if (lfirst) then - out_stat = param%out_stat - else - out_stat = 'APPEND' - end if - select case(out_stat) - case('APPEND') - open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) - case('NEW', 'REPLACE', 'UNKNOWN') - open(unit=LUN, file=param%discard_out, status=out_stat, form='FORMATTED', err=667, iomsg=errmsg) - case default - write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(out_stat)) - call util_exit(FAILURE) - end select - lfirst = .false. - if (param%lgr) call tp_discards%pv2v(param) - - write(LUN, HDRFMT) param%t, nsp, param%lbig_discard - do i = 1, nsp - write(LUN, NAMEFMT, err = 667, iomsg = errmsg) SUB, tp_discards%id(i), tp_discards%status(i) - write(LUN, VECFMT, err = 667, iomsg = errmsg) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i) - write(LUN, VECFMT, err = 667, iomsg = errmsg) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i) - end do - if (param%lbig_discard) then - if (param%lgr) then - allocate(pltemp, source = pl) - call pltemp%pv2v(param) - allocate(vh, source = pltemp%vh) - deallocate(pltemp) - else - allocate(vh, source = pl%vh) - end if - - write(LUN, NPLFMT) npl - do i = 1, npl - write(LUN, PLNAMEFMT, err = 667, iomsg = errmsg) pl%id(i), pl%Gmass(i), pl%radius(i) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i) - write(LUN, VECFMT, err = 667, iomsg = errmsg) vh(1, i), vh(2, i), vh(3, i) - end do - deallocate(vh) - end if - close(LUN) - end associate - - return - - 667 continue - write(*,*) "Error writing discard file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_write_discard - - - module subroutine io_write_frame_body(self, iu, param) - !! author: David A. Minton - !! - !! Write a frame of output of either test particle or massive body data to the binary output file - !! Note: If outputting to orbital elements, but sure that the conversion is done prior to calling this method - !! - !! Adapted from David E. Kaufmann's Swifter routine io_write_frame.f90 - !! Adapted from Hal Levison's Swift routine io_write_frame.F - implicit none - ! Arguments - class(swiftest_body), intent(in) :: self !! Swiftest particle 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 - ! Internals - character(len=STRMAX) :: errmsg - - 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%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) - write(iu, err = 667, iomsg = errmsg) self%xh(3, 1:n) - write(iu, err = 667, iomsg = errmsg) self%vh(1, 1:n) - write(iu, err = 667, iomsg = errmsg) self%vh(2, 1:n) - write(iu, err = 667, iomsg = errmsg) self%vh(3, 1:n) - end if - if ((param%out_form == EL) .or. (param%out_form == XVEL)) then - write(iu, err = 667, iomsg = errmsg) self%a(1:n) - write(iu, err = 667, iomsg = errmsg) self%e(1:n) - write(iu, err = 667, iomsg = errmsg) self%inc(1:n) * RAD2DEG - write(iu, err = 667, iomsg = errmsg) self%capom(1:n) * RAD2DEG - write(iu, err = 667, iomsg = errmsg) self%omega(1:n) * RAD2DEG - write(iu, err = 667, iomsg = errmsg) self%capm(1:n) * RAD2DEG - end if - select type(pl => self) - class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body - write(iu, err = 667, iomsg = errmsg) pl%Gmass(1:n) - if (param%lrhill_present) write(iu, err = 667, iomsg = errmsg) pl%rhill(1:n) - if (param%lclose) write(iu, err = 667, iomsg = errmsg) pl%radius(1:n) - if (param%lrotation) then - write(iu, err = 667, iomsg = errmsg) pl%Ip(1, 1:n) - write(iu, err = 667, iomsg = errmsg) pl%Ip(2, 1:n) - write(iu, err = 667, iomsg = errmsg) pl%Ip(3, 1:n) - write(iu, err = 667, iomsg = errmsg) pl%rot(1, 1:n) - write(iu, err = 667, iomsg = errmsg) pl%rot(2, 1:n) - write(iu, err = 667, iomsg = errmsg) pl%rot(3, 1:n) - end if - ! if (param%ltides) then - ! write(iu, err = 667, iomsg = errmsg) pl%k2(1:n) - ! write(iu, err = 667, iomsg = errmsg) pl%Q(1:n) - ! end if - end select - end associate - - return - 667 continue - write(*,*) "Error writing body frame: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_write_frame_body - - - module subroutine io_write_frame_cb(self, iu, param) - !! author: David A. Minton - !! - !! Write a frame of output of central body data to the binary output file - !! - !! Adapted from David E. Kaufmann's Swifter routine io_write_frame.f90 - !! Adapted from Hal Levison's Swift routine io_write_frame.F - implicit none - ! Arguments - class(swiftest_cb), intent(in) :: self !! Swiftest central body 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 - ! Internals - character(len=STRMAX) :: errmsg - - associate(cb => self) - write(iu, err = 667, iomsg = errmsg) cb%id - 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 - write(iu, err = 667, iomsg = errmsg) cb%j4rp4 - if (param%lrotation) then - write(iu, err = 667, iomsg = errmsg) cb%Ip(1) - write(iu, err = 667, iomsg = errmsg) cb%Ip(2) - write(iu, err = 667, iomsg = errmsg) cb%Ip(3) - write(iu, err = 667, iomsg = errmsg) cb%rot(1) - write(iu, err = 667, iomsg = errmsg) cb%rot(2) - write(iu, err = 667, iomsg = errmsg) cb%rot(3) - end if - ! if (param%ltides) then - ! write(iu, err = 667, iomsg = errmsg) cb%k2 - ! write(iu, err = 667, iomsg = errmsg) cb%Q - ! end if - end associate - - return - 667 continue - write(*,*) "Error writing central body frame: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_write_frame_cb - - module subroutine io_write_frame_system(self, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -2055,19 +1499,10 @@ module subroutine io_write_frame_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals logical, save :: lfirst = .true. !! Flag to determine if this is the first call of this method - class(swiftest_cb), allocatable :: cb !! Temporary local version of pl structure used for non-destructive conversions - class(swiftest_pl), allocatable :: pl !! Temporary local version of pl structure used for non-destructive conversions - class(swiftest_tp), allocatable :: tp !! Temporary local version of pl structure used for non-destructive conversions character(len=STRMAX) :: errmsg - integer(I4B) :: iu = BINUNIT !! Unit number for the output file to write frame to logical :: fileExists - allocate(cb, source = self%cb) - allocate(pl, source = self%pl) - allocate(tp, source = self%tp) - iu = BINUNIT - - param%nciu%id_chunk = pl%nbody + tp%nbody + param%nciu%id_chunk = self%pl%nbody + self%tp%nbody param%nciu%time_chunk = max(param%istep_dump / param%istep_out, 1) if (lfirst) then inquire(file=param%outfile, exist=fileExists) @@ -2091,12 +1526,7 @@ module subroutine io_write_frame_system(self, param) lfirst = .false. end if - ! Write out each data type frame - ! For NetCDF output, because we want to store the pseudovelocity separately from the true velocity, we need to do the orbital element conversion internally - call self%write_hdr(param%nciu, param) - call cb%write_frame(param%nciu, param) - call pl%write_frame(param%nciu, param) - call tp%write_frame(param%nciu, param) + call self%write_frame(param%nciu, param) return @@ -2105,37 +1535,4 @@ module subroutine io_write_frame_system(self, param) call util_exit(FAILURE) end subroutine io_write_frame_system - - module subroutine io_write_hdr_system(self, iu, param) ! t, npl, ntp, out_form, out_type) - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Write frame header to output binary file - !! - !! Adapted from David Adapted from David E. Kaufmann's Swifter routine io_write_hdr.f90 - !! Adapted from Hal Levison's Swift routine io_write_hdr.F - implicit none - ! Arguments - class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object - integer(I4B), intent(inout) :: iu !! Output file unit number - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - character(len=STRMAX) :: errmsg - - select case (param%out_type) - case (REAL4_TYPE) - write(iu, err = 667, iomsg = errmsg) real(param%t, kind=SP) - case (REAL8_TYPE) - write(iu, err = 667, iomsg = errmsg) param%t - end select - write(iu, err = 667, iomsg = errmsg) self%pl%nbody - write(iu, err = 667, iomsg = errmsg) self%tp%nbody - write(iu, err = 667, iomsg = errmsg) param%out_form - - return - - 667 continue - write(*,*) "Error writing header: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_write_hdr_system - end submodule s_io diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index c25566fe0..8f02b74f8 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -85,10 +85,8 @@ program swiftest_driver ioutput_t0 = int(t0 / dt / istep_out, kind=I8B) ioutput = ioutput_t0 ! Prevent duplicate frames from being written if this is a restarted run - if ((param%lrestart) .and. ((param%out_type == REAL8_TYPE) .or. param%out_type == REAL4_TYPE)) then - old_t_final = nbody_system%get_old_t_final_bin(param) - else if ((param%lrestart) .and. ((param%out_type == NETCDF_DOUBLE_TYPE) .or. param%out_type == NETCDF_FLOAT_TYPE)) then - old_t_final = nbody_system%get_old_t_final_netcdf(param) + if (param%lrestart) then + old_t_final = nbody_system%get_old_t_final(param) else old_t_final = t0 if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index bf209c477..b562def23 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -38,7 +38,6 @@ module encounter_classes procedure :: dealloc => encounter_util_dealloc_list !! Deallocates all allocatables procedure :: spill => encounter_util_spill_list !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) procedure :: resize => encounter_util_resize_list !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. - procedure :: write => encounter_io_write_list !! Write close encounter data to output binary file final :: encounter_util_final_list !! Finalize the encounter list - deallocates all allocatables end type encounter_list @@ -185,15 +184,6 @@ module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radi real(DP), dimension(:), intent(in) :: vh1, vh2 !! Swiftestcentric velocity vectors of the two encountering bodies end subroutine encounter_io_write_frame - module subroutine encounter_io_write_list(self, pl, encbody, param) - use swiftest_classes, only : swiftest_pl, swiftest_body, swiftest_parameters - implicit none - class(encounter_list), intent(in) :: self !! Swiftest encounter list object - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine encounter_io_write_list - module subroutine encounter_setup_aabb(self, n, n_last) implicit none class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 7fe65ffc9..5c3dc7adc 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -141,14 +141,6 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount logical :: lencounter !! Returns true if there is at least one close encounter end function rmvs_encounter_check_tp - module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2, enc_out) - implicit none - integer(I4B), intent(in) :: id1, id2 - real(DP), intent(in) :: t, Gmass1, Gmass2, radius1, radius2 - real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: enc_out - end subroutine rmvs_io_write_encounter - module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index f04346624..874bfdf31 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -46,12 +46,12 @@ module swiftest_classes character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles character(STRMAX) :: in_netcdf = NC_INFILE !! Name of system input file for NetCDF input - character(STRMAX) :: in_type = ASCII_TYPE !! Data representation type of input data files - character(STRMAX) :: in_form = XV !! Format of input data files (EL or XV) + character(STRMAX) :: in_type = "ASCII" !! Data representation type of input data files + character(STRMAX) :: in_form = "XV" !! Format of input data files ("EL" or "XV") integer(I4B) :: istep_out = -1 !! Number of time steps between binary outputs character(STRMAX) :: outfile = NETCDF_OUTFILE !! Name of output binary file - 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_type = "NETCDF_DOUBLE" !! 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 integer(I4B) :: istep_dump = -1 !! Number of time steps between dumps real(DP) :: rmin = -1.0_DP !! Minimum heliocentric radius for test particle @@ -61,14 +61,11 @@ module swiftest_classes character(STRMAX) :: qmin_coord = 'HELIO' !! Coordinate frame to use for qmin real(DP) :: qmin_alo = -1.0_DP !! Minimum semimajor axis for qmin real(DP) :: qmin_ahi = -1.0_DP !! Maximum semimajor axis for qmin - character(STRMAX) :: enc_out = "" !! Name of output file for encounters - character(STRMAX) :: discard_out = "" !! Name of output file for discards real(QP) :: MU2KG = -1.0_QP !! Converts mass units to grams real(QP) :: TU2S = -1.0_QP !! Converts time units to seconds 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) :: energy_out = "" !! Name of output energy and momentum report file character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" character(NAMELEN) :: encounter_check_plpl = "ADAPTIVE" !! Method used to compute pl-pl encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" character(NAMELEN) :: encounter_check_pltp = "ADAPTIVE" !! Method used to compute pl-tp encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" @@ -143,7 +140,6 @@ module swiftest_classes real(DP), dimension(NDIM) :: discard_vh !! The heliocentric velocity vector at the time of the particle's discard integer(I4B) :: discard_body_id !! The id of the other body involved in the discard (0 if no other body involved) contains - procedure :: read_in => io_read_in_particle_info !! Read in a particle information object from an open file procedure :: copy => util_copy_particle_info !! Copies one set of information object components into another, component-by-component procedure :: set_value => util_set_particle_info !! Sets one or more values of the particle information metadata object end type swiftest_particle_info @@ -155,12 +151,9 @@ module swiftest_classes !! An abstract superclass for a generic Swiftest object contains !! The minimal methods that all systems must have - procedure :: dump => io_dump_base !! Dump contents to file - procedure :: read_in => io_read_in_base !! Read in body initial conditions from a 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 - procedure :: write_particle_info_netcdf => netcdf_write_particle_info_base !! Dump contents of particle information metadata to file - generic :: write_frame => write_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments - generic :: write_particle_info => write_particle_info_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments + procedure :: read_in => io_read_in_base !! Read in body initial conditions from a file + procedure :: write_frame => 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 + procedure :: write_particle_info => netcdf_write_particle_info_base !! Dump contents of particle information metadata to file end type swiftest_base !******************************************************************************************************************************** @@ -193,10 +186,6 @@ module swiftest_classes 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_frame_bin => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body - procedure :: write_frame_bin => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body - generic :: write_frame => write_frame_bin !! Write a frame (either binary or NetCDF, using generic procedures) - generic :: read_frame => read_frame_bin !! Write a frame (either binary or NetCDF, using generic procedures) end type swiftest_cb !******************************************************************************************************************************** @@ -241,7 +230,6 @@ module swiftest_classes procedure :: v2pv => gr_vh2pv_body !! Converts from velocity to psudeovelocity for GR calculations using symplectic integrators procedure :: pv2v => gr_pv2vh_body !! Converts from psudeovelocity to velocity for GR calculations using symplectic integrators procedure :: read_frame_bin => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body - procedure :: write_frame_bin => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body procedure :: accel_obl => obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: el2xv => orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors procedure :: xv2el => orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements @@ -255,7 +243,6 @@ module swiftest_classes procedure :: sort => util_sort_body !! Sorts body arrays by a sortable componen procedure :: rearrange => util_sort_rearrange_body !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => util_spill_body !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - generic :: write_frame => write_frame_bin !! Add the generic write frame for Fortran binary files generic :: read_frame => read_frame_bin !! Add the generic read frame for Fortran binary files end type swiftest_body @@ -408,18 +395,14 @@ module swiftest_classes procedure :: compact_output => io_compact_output !! Prints out out terminal output when display_style is set to COMPACT 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_bin => 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 :: get_old_t_final_netcdf => netcdf_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. - procedure :: read_frame_bin => io_read_frame_system !! Read in a frame of input data from file - procedure :: write_frame_bin => io_write_frame_system !! Append a frame of output data to file - procedure :: read_frame_netcdf => netcdf_read_frame_system !! Read in a frame of input data from file + procedure :: get_old_t_final => netcdf_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. + procedure :: read_frame => netcdf_read_frame_system !! Read in a frame of input data from file procedure :: write_frame_netcdf => netcdf_write_frame_system !! Write a frame of input data from file - procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format - procedure :: read_hdr_netcdf => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format - procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format + procedure :: write_frame_system => io_write_frame_system !! Write a frame of input data from file + procedure :: read_hdr => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format + procedure :: write_hdr => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system - procedure :: read_particle_info_netcdf => netcdf_read_particle_info_system !! Read in particle metadata from file - procedure :: write_discard => io_write_discard !! Write out information about discarded test particles + procedure :: read_particle_info => netcdf_read_particle_info_system !! Read in particle metadata from file procedure :: obl_pot => obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body procedure :: finalize => setup_finalize_system !! Runs any finalization subroutines when ending the simulation. procedure :: initialize => setup_initialize_system !! Initialize the system from input files @@ -430,11 +413,7 @@ module swiftest_classes 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 - generic :: read_hdr => read_hdr_netcdf !! Generic method call for reading headers - generic :: read_frame => read_frame_bin, read_frame_netcdf !! Generic method call for reading a frame of output data - generic :: write_frame => write_frame_bin, write_frame_netcdf !! Generic method call for writing a frame of output data - generic :: read_particle_info => read_particle_info_netcdf !! Genereric method call for reading in the particle information metadata + generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data end type swiftest_nbody_system @@ -626,12 +605,6 @@ 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_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_base - module subroutine io_dump_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -645,13 +618,6 @@ module subroutine io_get_args(integrator, param_file_name, display_style) character(len=:), allocatable, intent(inout) :: display_style !! Style of the output display {"STANDARD", "COMPACT"}). Default is "STANDARD" end subroutine io_get_args - module function io_get_old_t_final_system(self, param) result(old_t_final) - implicit none - class(swiftest_nbody_system), intent(in) :: self - class(swiftest_parameters), intent(in) :: param - real(DP) :: old_t_final - end function io_get_old_t_final_system - module function io_get_token(buffer, ifirst, ilast, ierr) result(token) implicit none character(len=*), intent(in) :: buffer !! Input string buffer @@ -769,12 +735,6 @@ 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 subroutine io_read_in_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self @@ -789,18 +749,10 @@ module function io_read_frame_body(self, iu, param) result(ierr) integer(I4B) :: ierr !! Error code: returns 0 if the read is successful end function io_read_frame_body - module function io_read_frame_cb(self, iu, param) result(ierr) - implicit none - class(swiftest_cb), intent(inout) :: self !! Swiftest central body 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 - integer(I4B) :: ierr !! Error code: returns 0 if the read is successful - end function io_read_frame_cb - module function io_read_frame_system(self, iu, param) result(ierr) implicit none class(swiftest_nbody_system),intent(inout) :: self !! Swiftest system object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + integer(I4B), intent(inout) :: iu !! Unit number for the output file to read frame from class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters integer(I4B) :: ierr !! Error code: returns 0 if the read is successful end function io_read_frame_system @@ -816,39 +768,12 @@ module subroutine io_toupper(string) character(*), intent(inout) :: string !! String to make upper case end subroutine io_toupper - module subroutine io_write_discard(self, param) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_write_discard - - module subroutine io_write_frame_body(self, iu, param) - implicit none - class(swiftest_body), intent(in) :: self !! Swiftest body 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 io_write_frame_body - - module subroutine io_write_frame_cb(self, iu, param) - implicit none - class(swiftest_cb), intent(in) :: self !! Swiftest central body 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 io_write_frame_cb - module subroutine io_write_frame_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine io_write_frame_system - module subroutine io_write_hdr_system(self, iu, param) - implicit none - class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object - integer(I4B), intent(inout) :: iu !! Output file unit number - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine io_write_hdr_system - module subroutine kick_getacch_int_pl(self, param) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 49fa0b834..f8e5674e4 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -57,16 +57,6 @@ module swiftest_globals integer(I4B), parameter :: STRMAX = 512 !! Maximum size of character strings integer(I4B), parameter :: NAMELEN = 32 !! Maximum size of name strings - character(*), parameter :: ASCII_TYPE = 'ASCII' !! Symbolic name for ASCII file type - character(*), parameter :: REAL4_TYPE = 'REAL4' !! Symbolic name for binary file type REAL4 - character(*), parameter :: REAL8_TYPE = 'REAL8' !! Symbolic name for binary file type REAL8 - character(*), parameter :: NETCDF_FLOAT_TYPE = 'NETCDF_FLOAT' !! Symbolic name for binary file type REAL8 - character(*), parameter :: NETCDF_DOUBLE_TYPE = 'NETCDF_DOUBLE' !! Symbolic name for binary file type REAL8 - - character(*), parameter :: EL = 'EL' !! Symbolic name for binary output file contents for orbital elements - 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 = "Central Body" character(*), parameter :: PL_TYPE_NAME = "Massive Body" character(*), parameter :: TP_TYPE_NAME = "Test Particle" diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index aa599f94e..89824e4ec 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -202,9 +202,9 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_dim(self%ncid, TIME_DIMNAME, NF90_UNLIMITED, self%time_dimid), "netcdf_initialize_output nf90_def_dim time_dimid" ) ! 'y' dimension select case (param%out_type) - case(NETCDF_FLOAT_TYPE) + case("NETCDF_FLOAT") self%out_type = NF90_FLOAT - case(NETCDF_DOUBLE_TYPE) + case("NETCDF_DOUBLE") self%out_type = NF90_DOUBLE end select @@ -218,7 +218,7 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(self%ncid, PTYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%ptype_varid), "netcdf_initialize_output nf90_def_var ptype_varid" ) call check( nf90_def_var(self%ncid, STATUS_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%status_varid), "netcdf_initialize_output nf90_def_var status_varid" ) - if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_def_var(self%ncid, XHX_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%xhx_varid), "netcdf_initialize_output nf90_def_var xhx_varid" ) call check( nf90_def_var(self%ncid, XHY_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%xhy_varid), "netcdf_initialize_output nf90_def_var xhy_varid" ) call check( nf90_def_var(self%ncid, XHZ_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%xhz_varid), "netcdf_initialize_output nf90_def_var xhz_varid" ) @@ -237,7 +237,7 @@ module subroutine netcdf_initialize_output(self, param) end if - if ((param%out_form == EL) .or. (param%out_form == XVEL)) then + if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then call check( nf90_def_var(self%ncid, A_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%a_varid), "netcdf_initialize_output nf90_def_var a_varid" ) call check( nf90_def_var(self%ncid, E_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%e_varid), "netcdf_initialize_output nf90_def_var e_varid" ) call check( nf90_def_var(self%ncid, INC_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%inc_varid), "netcdf_initialize_output nf90_def_var inc_varid" ) @@ -381,7 +381,7 @@ module subroutine netcdf_open(self, param, readonly) call check( nf90_inq_varid(self%ncid, PTYPE_VARNAME, self%ptype_varid), "netcdf_open nf90_inq_varid ptype_varid" ) call check( nf90_inq_varid(self%ncid, GMASS_VARNAME, self%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" ) - if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_inq_varid(self%ncid, XHX_VARNAME, self%xhx_varid), "netcdf_open nf90_inq_varid xhx_varid" ) call check( nf90_inq_varid(self%ncid, XHY_VARNAME, self%xhy_varid), "netcdf_open nf90_inq_varid xhy_varid" ) call check( nf90_inq_varid(self%ncid, XHZ_VARNAME, self%xhz_varid), "netcdf_open nf90_inq_varid xhz_varid" ) @@ -408,7 +408,7 @@ module subroutine netcdf_open(self, param, readonly) end if end if - if ((param%out_form == EL) .or. (param%out_form == XVEL)) then + if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then call check( nf90_inq_varid(self%ncid, A_VARNAME, self%a_varid), "netcdf_open nf90_inq_varid a_varid" ) call check( nf90_inq_varid(self%ncid, E_VARNAME, self%e_varid), "netcdf_open nf90_inq_varid e_varid" ) call check( nf90_inq_varid(self%ncid, INC_VARNAME, self%inc_varid), "netcdf_open nf90_inq_varid inc_varid" ) @@ -531,7 +531,7 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) call check( nf90_inquire_dimension(iu%ncid, iu%str_dimid, len=str_max), "netcdf_read_frame_system nf90_inquire_dimension str_dimid" ) ! First filter out only the id slots that contain valid bodies - if (param%in_form == XV) then + if (param%in_form == "XV") then call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system filter pass nf90_getvar xhx_varid" ) else call check( nf90_get_var(iu%ncid, iu%a_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system filter pass nf90_getvar a_varid" ) @@ -572,7 +572,7 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) end select ! Now read in each variable and split the outputs by body type - if ((param%in_form == XV) .or. (param%in_form == XVEL)) then + if ((param%in_form == "XV") .or. (param%in_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar xhx_varid" ) if (npl > 0) pl%xh(1,:) = pack(rtemp, plmask) if (ntp > 0) tp%xh(1,:) = pack(rtemp, tpmask) @@ -612,7 +612,7 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) end if end if - if ((param%in_form == EL) .or. (param%in_form == XVEL)) then + if ((param%in_form == "EL") .or. (param%in_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%a_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar a_varid" ) if (.not.allocated(pl%a)) allocate(pl%a(npl)) if (.not.allocated(tp%a)) allocate(tp%a(ntp)) @@ -995,7 +995,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_XHX_VARNAME, iu%origin_xhx_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_xhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar origin_xhx_varid" ) - else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar xhx_varid" ) else rtemp_arr(1,:) = 0._DP @@ -1004,7 +1004,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_XHY_VARNAME, iu%origin_xhy_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_xhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar origin_xhy_varid" ) - else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%xhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar xhx_varid" ) else rtemp_arr(2,:) = 0._DP @@ -1013,7 +1013,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_XHZ_VARNAME, iu%origin_xhz_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_xhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar origin_xhz_varid" ) - else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%xhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar xhz_varid" ) else rtemp_arr(3,:) = 0._DP @@ -1029,7 +1029,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_VHX_VARNAME, iu%origin_vhx_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_vhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar origin_vhx_varid" ) - else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar vhx_varid" ) else rtemp_arr(1,:) = 0._DP @@ -1038,7 +1038,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_VHY_VARNAME, iu%origin_vhy_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_vhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar origin_vhy_varid" ) - else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar vhy_varid" ) else rtemp_arr(2,:) = 0._DP @@ -1047,7 +1047,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_VHZ_VARNAME, iu%origin_vhz_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_vhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar origin_vhz_varid" ) - else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar vhz_varid" ) else rtemp_arr(3,:) = 0._DP @@ -1202,7 +1202,7 @@ module subroutine netcdf_write_frame_base(self, iu, param) !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity if (param%lgr) call gr_pseudovel2vel(param, self%mu(j), self%xh(:, j), self%vh(:, j), vh(:)) - if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then call check( nf90_put_var(iu%ncid, iu%xhx_varid, self%xh(1, j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var xhx_varid" ) call check( nf90_put_var(iu%ncid, iu%xhy_varid, self%xh(2, j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var xhy_varid" ) call check( nf90_put_var(iu%ncid, iu%xhz_varid, self%xh(3, j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var xhz_varid" ) @@ -1221,7 +1221,7 @@ module subroutine netcdf_write_frame_base(self, iu, param) end if end if - if ((param%out_form == EL) .or. (param%out_form == XVEL)) then + if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then if (param%lgr) then !! For GR-enabled runs, use the true value of velocity computed above call orbel_xv2el(self%mu(j), self%xh(1,j), self%xh(2,j), self%xh(3,j), & vh(1), vh(2), vh(3), & diff --git a/src/rmvs/rmvs_io.f90 b/src/rmvs/rmvs_io.f90 deleted file mode 100644 index 4d04dc150..000000000 --- a/src/rmvs/rmvs_io.f90 +++ /dev/null @@ -1,56 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (rmvs_classes) s_rmvs_io - use swiftest -contains - - module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & - xh1, xh2, vh1, vh2, enc_out) - !! author: David A. Minton - !! - !! Write close encounter data from RMVS to output binary files - !! There is no direct file output from this subroutine - !! - !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 - !! Adapted from Hal Levison's Swift routine io_write_encounter.f - implicit none - ! Arguments - integer(I4B), intent(in) :: id1, id2 - real(DP), intent(in) :: t, Gmass1, Gmass2, radius1, radius2 - real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: enc_out - ! Internals - logical , save :: lfirst = .true. - integer(I4B) :: ierr - - if (enc_out == "") return - - open(unit = LUN, file = enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) - if ((ierr /= 0) .and. lfirst) then - open(unit = LUN, file = enc_out, status = 'NEW', form = 'UNFORMATTED', iostat = ierr) - end if - if (ierr /= 0) then - write(*, *) "Swiftest Error:" - write(*, *) " Unable to open binary encounter file" - call util_exit(FAILURE) - end if - lfirst = .false. - call encounter_io_write_frame(LUN, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) - close(unit = LUN, iostat = ierr) - if (ierr /= 0) then - write(*, *) "Swiftest Error:" - write(*, *) " Unable to close binary encounter file" - call util_exit(FAILURE) - end if - - return - end subroutine rmvs_io_write_encounter - -end submodule s_rmvs_io \ No newline at end of file diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 25285899f..7c39614e1 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -561,17 +561,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) a, peri, capm, tperi) r2 = dot_product(xpc(:, i), xpc(:, i)) if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) - if (param%enc_out /= "") then - id1 = pl%id(ipleP) - rpl = pl%radius(ipleP) - xh1(:) = pl%inner(inner_index)%x(:, ipleP) - vh1(:) = pl%inner(inner_index)%v(:, ipleP) - id2 = tp%id(i) - xh2(:) = xpc(:, i) + xh1(:) - vh2(:) = xpc(:, i) + vh1(:) - call rmvs_io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, & - xh1(:), xh2(:), vh1(:), vh2(:), param%enc_out) - end if + ! TODO: write NetCDF encounter output writer if (tp%lperi(i)) then if (peri < tp%peri(i)) then tp%peri(i) = peri diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 3d0943d95..ae157eaae 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -91,7 +91,7 @@ module subroutine setup_finalize_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters associate(system => self) - if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + if ((param%out_type == "NETCDF_FLOAT") .or. (param%out_type == "NETCDF_DOUBLE")) then call param%nciu%close() end if end associate @@ -148,7 +148,7 @@ module subroutine setup_initialize_system(self, param) call system%set_msys() call pl%set_mu(cb) call tp%set_mu(cb) - if (param%in_form == EL) then + if (param%in_form == "EL") then call pl%el2xv(cb) call tp%el2xv(cb) end if diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f60b91a28..f00b222cb 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -361,7 +361,7 @@ module subroutine symba_discard_pl(self, system, param) associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list) call pl%vb2vh(system%cb) call pl%xh2xb(system%cb) - call plplenc_list%write(pl, pl, param) + !call plplenc_list%write(pl, pl, param) TODO: write the encounter list writer for NetCDF call symba_discard_nonplpl(self, system, param) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index d9a48b52a..d5dd06308 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -174,85 +174,20 @@ module subroutine symba_io_write_discard(self, param) class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: iadd, isub, j, nsub, nadd - logical, save :: lfirst = .true. - character(*), parameter :: HDRFMT = '(E23.16, 1X, I8, 1X, L1)' - character(*), parameter :: NAMEFMT = '(A, 2(1X, I8))' - character(*), parameter :: VECFMT = '(3(E23.16, 1X))' - character(*), parameter :: NPLFMT = '(I8)' - character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' - character(STRMAX) :: errmsg, out_stat associate(pl => self%pl, npl => self%pl%nbody, pl_adds => self%pl_adds) - if (self%tp_discards%nbody > 0) call io_write_discard(self, param) + if (self%tp_discards%nbody > 0) call self%tp_discards%write_particle_info(param%nciu, param) select type(pl_discards => self%pl_discards) class is (symba_merger) if (pl_discards%nbody == 0) return - ! Record the discarded body metadata information to file - if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - call pl_discards%write_particle_info(param%nciu, param) - end if - - if (param%discard_out == "") return - if (lfirst) then - out_stat = param%out_stat - else - out_stat = 'APPEND' - end if - select case(out_stat) - case('APPEND') - open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) - case('NEW', 'REPLACE', 'UNKNOWN') - open(unit=LUN, file=param%discard_out, status=param%out_stat, form='FORMATTED', err=667, iomsg=errmsg) - case default - write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select - lfirst = .false. - if (param%lgr) then - call pl_discards%pv2v(param) - call pl_adds%pv2v(param) - end if - - write(LUN, HDRFMT, err=667, iomsg=errmsg) param%t, pl_discards%nbody, param%lbig_discard - iadd = 1 - isub = 1 - do while (iadd <= pl_adds%nbody) - nadd = pl_adds%ncomp(iadd) - nsub = pl_discards%ncomp(isub) - do j = 1, nadd - if (iadd <= pl_adds%nbody) then - write(LUN, NAMEFMT, err=667, iomsg=errmsg) ADD, pl_adds%id(iadd), pl_adds%status(iadd) - write(LUN, VECFMT, err=667, iomsg=errmsg) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) - write(LUN, VECFMT, err=667, iomsg=errmsg) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) - else - exit - end if - iadd = iadd + 1 - end do - do j = 1, nsub - if (isub <= pl_discards%nbody) then - write(LUN,NAMEFMT,err=667,iomsg=errmsg) SUB, pl_discards%id(isub), pl_discards%status(isub) - write(LUN,VECFMT,err=667,iomsg=errmsg) pl_discards%xh(1,isub), pl_discards%xh(2,isub), pl_discards%xh(3,isub) - write(LUN,VECFMT,err=667,iomsg=errmsg) pl_discards%vh(1,isub), pl_discards%vh(2,isub), pl_discards%vh(3,isub) - else - exit - end if - isub = isub + 1 - end do - end do - - close(LUN) + call pl_discards%write_particle_info(param%nciu, param) end select end associate return - 667 continue - write(*,*) "Error writing discard file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) end subroutine symba_io_write_discard end submodule s_symba_io diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 8196a5f77..a9755d0d4 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -89,7 +89,7 @@ module subroutine whm_setup_initialize_system(self, param) call self%tp_discards%setup(0, param) call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) - if (param%lgr .and. ((param%in_type == REAL8_TYPE) .or. (param%in_type == REAL4_TYPE))) then !! pseudovelocity conversion for NetCDF input files is handled by NetCDF routines + if (param%lgr .and. param%in_type == "ASCII") then !! pseudovelocity conversion for NetCDF input files is handled by NetCDF routines call self%pl%v2pv(param) call self%tp%v2pv(param) end if