diff --git a/.gitignore b/.gitignore index 67bc9fa6b..fe0a24d5d 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ !README.md !README.swifter !*.in +dump* !example/cleanup !example/swifter_symba_omp !Makefile diff --git a/src/io/io.f90 b/src/io/io.f90 index 4aa3b88a3..e2171b169 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -157,8 +157,12 @@ module subroutine io_dump_swiftest(self, param) 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 + end select call self%write_frame(iu, param) - close(LUN, err = 667, iomsg = errmsg) + close(iu, err = 667, iomsg = errmsg) return @@ -204,13 +208,15 @@ module subroutine io_dump_system(self, param) dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx))) dump_param%out_form = XV + dump_param%in_form = XV dump_param%out_stat = 'APPEND' + dump_param%in_type = REAL8_TYPE dump_param%T0 = param%t call dump_param%dump(param_file_name) call self%cb%dump(dump_param) - if (self%pl%nbody > 0) call self%pl%dump(dump_param) - if (self%tp%nbody > 0) call self%tp%dump(dump_param) + call self%pl%dump(dump_param) + call self%tp%dump(dump_param) idx = idx + 1 if (idx > NDUMPFILES) idx = 1 @@ -293,6 +299,45 @@ module function io_get_args(integrator, param_file_name) result(ierr) end function 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), parameter :: LUN = 76 + 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 + 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 !! @@ -404,6 +449,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("IN_TYPE") call io_toupper(param_value) param%in_type = param_value + case ("IN_FORM") + call io_toupper(param_value) + param%in_form = param_value case ("ISTEP_OUT") read(param_value, *) param%istep_out case ("BIN_OUT") @@ -606,6 +654,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "BIN_OUT = ",trim(adjustl(param%outfile)) write(*,*) "OUT_TYPE = ",trim(adjustl(param%out_type)) write(*,*) "OUT_FORM = ",trim(adjustl(param%out_form)) + write(*,*) "IN_FORM = ",trim(adjustl(param%in_form)) write(*,*) "OUT_STAT = ",trim(adjustl(param%out_stat)) write(*,*) "ISTEP_DUMP = ",param%istep_dump write(*,*) "CHK_CLOSE = ",param%lclose @@ -719,9 +768,11 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "T0"; write(param_value,Rfmt) param%t0; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "TSTOP"; write(param_value, Rfmt) param%tstop; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "DT"; write(param_value, Rfmt) param%dt; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "CB_IN"; write(param_value, Afmt) trim(adjustl(param%incbfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "PL_IN"; write(param_value, Afmt) trim(adjustl(param%inplfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "TP_in"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "TP_IN"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "IN_TYPE"; write(param_value, Afmt) trim(adjustl(param%in_type)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "IN_FORM"; write(param_value, Afmt) trim(adjustl(param%in_form)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) if (param%istep_out > 0) then write(param_name, Afmt) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "BIN_OUT"; write(param_value, Afmt) trim(adjustl(param%outfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) @@ -749,9 +800,11 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "DU2M"; write(param_value, Rfmt) param%DU2M; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "RHILL_PRESENT"; write(param_value, Lfmt) param%lrhill_present; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "EXTRA_FORCE"; write(param_value, Lfmt) param%lextra_force; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "BIG_DISCARD"; write(param_value, Lfmt) param%lbig_discard; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "DISCARD_OUT"; write(param_value, Afmt) trim(adjustl(param%discard_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + if (param%discard_out /= "") write(param_name, Afmt) "BIG_DISCARD"; write(param_value, Lfmt) param%lbig_discard; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "CHK_CLOSE"; write(param_value, Lfmt) param%lclose; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "ENERGY"; write(param_value, Lfmt) param%lenergy; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + if (param%lenergy) write(param_name, Afmt) "ENERGY_OUT"; write(param_value, Afmt) trim(adjustl(param%energy_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "GR"; write(param_value, Lfmt) param%lgr; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) @@ -781,7 +834,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) end subroutine io_param_writer - module subroutine io_read_body_in(self, param) + module subroutine io_read_in_body(self, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Read in either test particle or massive body data @@ -795,12 +848,12 @@ module subroutine io_read_body_in(self, param) ! Internals integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: iu = LUN - integer(I4B) :: i, nbody + integer(I4B) :: nbody logical :: is_ascii, is_pl character(len=:), allocatable :: infile real(DP) :: t - real(QP) :: val character(STRMAX) :: errmsg + integer(I4B) :: ierr ! Select the appropriate polymorphic class (test particle or massive body) select type(self) @@ -817,63 +870,32 @@ module subroutine io_read_body_in(self, param) case(ASCII_TYPE) 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) - if (nbody > 0) then - do i = 1, nbody - select type(self) - class is (swiftest_pl) - if (param%lrhill_present) then - read(iu, *, err = 667, iomsg = errmsg) self%id(i), val, self%rhill(i) - else - read(iu, *, err = 667, iomsg = errmsg) self%id(i), val - end if - self%Gmass(i) = real(val, kind=DP) - self%mass(i) = real(val / param%GU, kind=DP) - if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i) - class is (swiftest_tp) - read(iu, *, err = 667, iomsg = errmsg) self%id(i) - end select - 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) - select type (self) - class is (swiftest_pl) - if (param%lrotation) then - read(iu, *, err = 667, iomsg = errmsg) self%Ip(1, i), self%Ip(2, i), self%Ip(3, i) - read(iu, *, err = 667, iomsg = errmsg) self%rot(1, i), self%rot(2, i), self%rot(3, i) - end if - if (param%ltides) then - read(iu, *, err = 667, iomsg = errmsg) self%k2(i) - read(iu, *, err = 667, iomsg = errmsg) self%Q(i) - end if - end select - self%status(i) = ACTIVE - self%lmask(i) = .true. - end do - end if - case (REAL4_TYPE, REAL8_TYPE) !, SWIFTER_REAL4_TYPE, SWIFTER_REAL8_TYPE) + 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 - call self%setup(nbody, param) - if (nbody > 0) then - call self%read_frame(iu, param, XV) - self%status(:) = ACTIVE - self%lmask(:) = .true. - end if case default write(errmsg,*) trim(adjustl(param%in_type)) // ' is an unrecognized file type' goto 667 end select + + call self%setup(nbody, param) + ierr = 0 + if (nbody > 0) then + ierr = self%read_frame(iu, param) + self%status(:) = ACTIVE + self%lmask(:) = .true. + end if close(iu, err = 667, iomsg = errmsg) - return + if (ierr == 0) return 667 continue write(*,*) 'Error reading in initial conditions file: ',trim(adjustl(errmsg)) return - end subroutine io_read_body_in + end subroutine io_read_in_body - module subroutine io_read_cb_in(self, param) + module subroutine io_read_in_cb(self, param) !! author: David A. Minton !! !! Reads in central body data @@ -888,6 +910,7 @@ module subroutine io_read_cb_in(self, param) integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: iu = LUN character(len=STRMAX) :: errmsg + integer(I4B) :: ierr if (param%in_type == 'ASCII') then open(unit = iu, file = param%incbfile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg) @@ -903,19 +926,28 @@ module subroutine io_read_cb_in(self, param) end if else open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - call self%read_frame(iu, param, XV) + ierr = self%read_frame(iu, param) end if close(iu, err = 667, iomsg = errmsg) - if (self%j2rp2 /= 0.0_DP) param%loblatecb = .true. - if (param%rmin < 0.0) param%rmin = self%radius - + if (ierr == 0) then + + if (self%j2rp2 /= 0.0_DP) param%loblatecb = .true. + if (param%rmin < 0.0) param%rmin = self%radius + + select type(cb => self) + class is (symba_cb) + cb%M0 = cb%mass + cb%R0 = cb%radius + cb%L0(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) + end select + end if return 667 continue write(*,*) "Error reading central body file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine io_read_cb_in + end subroutine io_read_in_cb function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & @@ -969,7 +1001,7 @@ function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & end function io_read_encounter - module subroutine io_read_frame_body(self, iu, param, form) + module function io_read_frame_body(self, iu, param) result(ierr) !! author: David A. Minton !! !! Reads a frame of output of either test particle or massive body data from a binary output file @@ -981,65 +1013,131 @@ module subroutine io_read_frame_body(self, iu, param, form) class(swiftest_body), intent(inout) :: self !! Swiftest particle 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") + ! Result + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals character(len=STRMAX) :: errmsg + integer(I4B) :: i + real(QP) :: val + + if (self%nbody == 0) return + + 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) - read(iu, err = 667, iomsg = errmsg) self%id(:) - !read(iu, err = 667, iomsg = errmsg) self%name(1:n) - select case (form) - case (EL) + + 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)) if (.not.allocated(self%capom)) allocate(self%capom(n)) if (.not.allocated(self%omega)) allocate(self%omega(n)) if (.not.allocated(self%capm)) allocate(self%capm(n)) - 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(:) - 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, :) - 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%rot(1, :) - read(iu, err = 667, iomsg = errmsg) pl%rot(2, :) - read(iu, err = 667, iomsg = errmsg) pl%rot(3, :) - 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, :) - end if - if (param%ltides) then - read(iu, err = 667, iomsg = errmsg) pl%k2(1:n) - read(iu, err = 667, iomsg = errmsg) pl%Q(1:n) - end if + end if + + select case(param%in_type) + case (REAL4_TYPE, REAL8_TYPE) + read(iu, err = 667, iomsg = errmsg) self%id(:) + + 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(1:n) + read(iu, err = 667, iomsg = errmsg) pl%Q(1:n) + end if + end select + + case (ASCII_TYPE) + do i = 1, n + + select type(self) + class is (swiftest_pl) + if (param%lrhill_present) then + read(iu, *, err = 667, iomsg = errmsg) self%id(i), val, self%rhill(i) + else + read(iu, *, err = 667, iomsg = errmsg) self%id(i), val + end if + self%Gmass(i) = real(val, kind=DP) + self%mass(i) = real(val / param%GU, kind=DP) + if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i) + class is (swiftest_tp) + read(iu, *, err = 667, iomsg = errmsg) self%id(i) + end select + + select case(param%in_form) + 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) + 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 + + select type (self) + class is (swiftest_pl) + if (param%lrotation) then + read(iu, *, err = 667, iomsg = errmsg) self%Ip(1, i), self%Ip(2, i), self%Ip(3, i) + read(iu, *, err = 667, iomsg = errmsg) self%rot(1, i), self%rot(2, i), self%rot(3, i) + end if + if (param%ltides) then + read(iu, *, err = 667, iomsg = errmsg) self%k2(i) + read(iu, *, err = 667, iomsg = errmsg) self%Q(i) + end if + end select + + end do end select end associate + ierr = 0 return 667 continue - write(*,*) "Error reading central body file: " // trim(adjustl(errmsg)) + select type (self) + class is (swiftest_pl) + write(*,*) "Error reading massive body file: " // trim(adjustl(errmsg)) + class is (swiftest_tp) + write(*,*) "Error reading test particle file: " // trim(adjustl(errmsg)) + class default + write(*,*) "Error reading body file: " // trim(adjustl(errmsg)) + end select call util_exit(FAILURE) - end subroutine io_read_frame_body + end function io_read_frame_body - module subroutine io_read_frame_cb(self, iu, param, form) + 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 @@ -1051,69 +1149,87 @@ module subroutine io_read_frame_cb(self, iu, param, form) 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 - character(*), intent(in) :: form !! Input format code ("XV" or "EL") + ! Result + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals character(len=STRMAX) :: errmsg - read(iu, err = 667, iomsg = errmsg) self%id !read(iu, err = 667, iomsg = errmsg) self%name + read(iu, err = 667, iomsg = errmsg) self%id 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(:) - read(iu, err = 667, iomsg = errmsg) self%rot(:) + 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 file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine io_read_frame_cb + end function io_read_frame_cb - module subroutine io_read_frame_system(self, iu, param, form) + 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_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 - character(*), intent(in) :: form !! Input format code ("XV" or "EL") + ! Result + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful ! Internals - logical, save :: lfirst = .true. character(len=STRMAX) :: errmsg - integer(I4B) :: ierr - iu = BINUNIT - if (lfirst) then - open(unit = iu, file = param%outfile, status = 'OLD', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - lfirst = .false. - end if ierr = io_read_hdr(iu, param%t, self%pl%nbody, self%tp%nbody, param%out_form, param%out_type) + if (is_iostat_end(ierr)) return ! Reached the end of the frames + if (ierr /= 0) then - write(*, *) "Swiftest error:" - write(*, *) " Binary output file already exists or cannot be accessed" - return + 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 + + if (param%in_form == EL) then + call self%pl%el2xv(self%cb) + call self%tp%el2xv(self%cb) end if - call self%cb%read_frame(iu, param, form) - call self%pl%read_frame(iu, param, form) - call self%tp%read_frame(iu, param, form) return 667 continue write(*,*) "Error reading system frame: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_read_frame_system + end function io_read_frame_system function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) @@ -1127,7 +1243,7 @@ function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) ! Arguments integer(I4B), intent(in) :: iu integer(I4B), intent(out) :: npl, ntp - character(*), intent(out) :: out_form + character(*), intent(out) :: out_form real(DP), intent(out) :: t character(*), intent(in) :: out_type ! Result @@ -1137,29 +1253,30 @@ function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) character(len=STRMAX) :: errmsg select case (out_type) - case (REAL4_TYPE, SWIFTER_REAL4_TYPE) - read(iu, iostat = ierr, err = 667, iomsg = errmsg) ttmp, npl, ntp, out_form - if (ierr /= 0) return + case (REAL4_TYPE) + read(iu, iostat = ierr, err = 667, iomsg = errmsg, end = 333) ttmp t = ttmp - case (REAL8_TYPE, SWIFTER_REAL8_TYPE) - read(iu, iostat = ierr, err = 667, iomsg = errmsg) t - 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 + 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_param_in(self, param_file_name) + module subroutine io_read_in_param(self, param_file_name) !! author: David A. Minton !! !! Read in parameters for the integration @@ -1191,7 +1308,7 @@ module subroutine io_read_param_in(self, param_file_name) 667 continue write(*,*) "Error reading parameter file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine io_read_param_in + end subroutine io_read_in_param module subroutine io_toupper(string) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 5e28452e0..3a274e596 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -19,6 +19,7 @@ program swiftest_driver integer(I8B) :: iout !! Output cadence counter integer(I8B) :: nloops !! Number of steps to take in the simulation integer(I4B) :: iu !! Unit number of binary file + real(DP) :: old_t_final = 0.0_DP !! Output time at which writing should start, in order to prevent duplicate lines being written for restarts ierr = io_get_args(integrator, param_file_name) if (ierr /= 0) then @@ -35,7 +36,7 @@ program swiftest_driver param%integrator = integrator call setup_construct_system(nbody_system, param) - call param%read_from_file(param_file_name) + call param%read_in(param_file_name) associate(t => param%t, & t0 => param%t0, & @@ -50,8 +51,15 @@ program swiftest_driver iout = istep_out idump = istep_dump nloops = ceiling(tstop / dt, kind=I8B) - if (istep_out > 0) call nbody_system%write_frame(iu, param) - call nbody_system%dump(param) + ! Prevent duplicate frames from being written if this is a restarted run + if (param%lrestart) then + old_t_final = nbody_system%get_old_t_final(param) + else + old_t_final = t0 + if (istep_out > 0) call nbody_system%write_frame(iu, param) + call nbody_system%dump(param) + end if + !> Define the maximum number of threads nthreads = 1 ! In the *serial* case @@ -73,7 +81,7 @@ program swiftest_driver if (istep_out > 0) then iout = iout - 1 if (iout == 0) then - call nbody_system%write_frame(iu, param) + if (t > old_t_final) call nbody_system%write_frame(iu, param) iout = istep_out end if end if diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index aed6f23b0..9a674c830 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -24,7 +24,8 @@ module swiftest_classes character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body 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_type = ASCII_TYPE !! Format of input data files + 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) integer(I4B) :: istep_out = -1 !! Number of time steps between binary outputs character(STRMAX) :: outfile = BIN_OUTFILE !! Name of output binary file character(STRMAX) :: out_type = REAL8_TYPE !! Binary format of output file @@ -78,10 +79,10 @@ module swiftest_classes logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect logical :: lyorp = .false. !! Turn on YORP effect contains - procedure :: reader => io_param_reader - procedure :: writer => io_param_writer - procedure :: dump => io_dump_param - procedure :: read_from_file => io_read_param_in + procedure :: reader => io_param_reader + procedure :: writer => io_param_writer + procedure :: dump => io_dump_param + procedure :: read_in => io_read_in_param end type swiftest_parameters !******************************************************************************************************************************** @@ -93,7 +94,6 @@ module swiftest_classes contains !! The minimal methods that all systems must have procedure :: dump => io_dump_swiftest - procedure(abstract_initialize), deferred :: initialize procedure(abstract_read_frame), deferred :: read_frame procedure(abstract_write_frame), deferred :: write_frame end type swiftest_base @@ -128,7 +128,7 @@ 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 :: initialize => io_read_cb_in !! I/O routine for reading in central body data + procedure :: read_in => io_read_in_cb !! I/O routine for reading in central body data procedure :: read_frame => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body procedure :: write_frame => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body end type swiftest_cb @@ -174,7 +174,7 @@ module swiftest_classes procedure :: drift => drift_body !! Loop through bodies and call Danby drift routine on heliocentric variables 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 :: initialize => io_read_body_in !! Read in body initial conditions from a file + procedure :: read_in => io_read_in_body !! Read in body initial conditions from a file procedure :: read_frame => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body procedure :: write_frame => 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 @@ -306,6 +306,7 @@ module swiftest_classes procedure :: discard => discard_system !! Perform a discard step on the system procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time procedure :: dump => io_dump_system !! Dump the state of the system to a file + procedure :: get_old_t_final => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file procedure :: write_discard => io_write_discard !! Write out information about discarded test particles procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file @@ -356,12 +357,6 @@ subroutine abstract_accel(self, system, param, t, lbeg) logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step end subroutine abstract_accel - subroutine abstract_initialize(self, param) - import swiftest_base, swiftest_parameters - class(swiftest_base), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine abstract_initialize - subroutine abstract_kick_body(self, system, param, t, dt, lbeg) import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP implicit none @@ -373,13 +368,13 @@ subroutine abstract_kick_body(self, system, param, t, dt, lbeg) logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine abstract_kick_body - subroutine abstract_read_frame(self, iu, param, form) + function abstract_read_frame(self, iu, param) result(ierr) import DP, I4B, swiftest_base, swiftest_parameters class(swiftest_base), intent(inout) :: self !! Swiftest base 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") - end subroutine abstract_read_frame + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + end function abstract_read_frame subroutine abstract_set_mu(self, cb) import swiftest_body, swiftest_cb @@ -579,6 +574,13 @@ module function io_get_args(integrator, param_file_name) result(ierr) integer(I4B) :: ierr !! I/O error code end function 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 @@ -610,47 +612,47 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine io_param_writer - module subroutine io_read_body_in(self, param) + module subroutine io_read_in_body(self, param) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_read_body_in + end subroutine io_read_in_body - module subroutine io_read_cb_in(self, param) + module subroutine io_read_in_cb(self, param) implicit none class(swiftest_cb), intent(inout) :: self !! Swiftest central body object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_read_cb_in + end subroutine io_read_in_cb - module subroutine io_read_param_in(self, param_file_name) + module subroutine io_read_in_param(self, param_file_name) implicit none class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) - end subroutine io_read_param_in + end subroutine io_read_in_param - module subroutine io_read_frame_body(self, iu, param, form) + module function io_read_frame_body(self, iu, param) result(ierr) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest 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 - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - end subroutine io_read_frame_body + class(swiftest_body), intent(inout) :: self !! Swiftest 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_body - module subroutine io_read_frame_cb(self, iu, param, form) + 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 - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - end subroutine io_read_frame_cb + 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 subroutine io_read_frame_system(self, iu, param, form) + 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 class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - end subroutine io_read_frame_system + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + end function io_read_frame_system module subroutine io_write_discard(self, param) implicit none diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 454f454f5..17c5d9b19 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -107,7 +107,7 @@ module swiftest_globals character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] - character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.dat', 'dump_param2.dat'] + character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.in', 'dump_param2.in'] !> Default file names that can be changed by the user in the parameters file character(*), parameter :: CB_INFILE = 'cb.in' diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ea15bf1fe..7ee72fd3b 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -133,16 +133,18 @@ module subroutine setup_initialize_system(self, param) ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - call self%cb%initialize(param) - call self%pl%initialize(param) - call self%tp%initialize(param) + + call self%cb%read_in(param) + call self%pl%read_in(param) + call self%tp%read_in(param) call self%validate_ids(param) call self%set_msys() call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) call self%pl%eucl_index() if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) + self%pl%lfirst = param%lfirstkick + self%tp%lfirst = param%lfirstkick return end subroutine setup_initialize_system diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index bf70f15dd..2bf7b2a90 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -92,7 +92,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) system%Mescape = system%Mescape + pl%mass(ipl) do i = 1, pl%nbody if (i == ipl) cycle - pe = pe - pl%mass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i)) + pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i)) end do Ltot(:) = 0.0_DP diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 5ab525e23..06472aae9 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -276,7 +276,7 @@ module subroutine symba_io_read_particle(system, param) associate(npl => pl%nbody, ntp => tp%nbody) do lmatch = .false. - read(LUN, err = 667, iomsg = errmsg) id + read(LUN, err = 667, iomsg = errmsg, end = 333) id if (idx == cb%id) then read(LUN, err = 667, iomsg = errmsg) cb%info @@ -308,6 +308,7 @@ module subroutine symba_io_read_particle(system, param) end select end select + 333 continue return 667 continue diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index ed5985c0b..fe2d74b1a 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -28,10 +28,8 @@ module subroutine symba_step_system(self, param, t, dt) call self%reset(param) lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt, 0) if (lencounter) then - tp%lfirst = pl%lfirst call self%interp(param, t, dt) - pl%lfirst = .true. - tp%lfirst = .true. + param%lfirstkick = .true. else self%irec = -1 call helio_step_system(self, param, t, dt) @@ -266,6 +264,9 @@ module subroutine symba_step_reset_system(self, param) call system%pl_adds%setup(0, param) call system%pl_discards%setup(0, param) + + tp%lfirst = param%lfirstkick + pl%lfirst = param%lfirstkick end associate end select end select diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index a80c8e646..08b5728e8 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -455,15 +455,15 @@ module subroutine symba_util_rearray_pl(self, system, param) end if ! Reset all of the status flags for this body - where(pl%status(:) /= INACTIVE) - pl%status(:) = ACTIVE - pl%ldiscard(:) = .false. - pl%lcollision(:) = .false. - pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY - pl%lmask(:) = .true. + where(pl%status(1:pl%nbody) /= INACTIVE) + pl%status(1:pl%nbody) = ACTIVE + pl%ldiscard(1:pl%nbody) = .false. + pl%lcollision(1:pl%nbody) = .false. + pl%lmtiny(1:pl%nbody) = pl%Gmass(1:pl%nbody) > param%GMTINY + pl%lmask(1:pl%nbody) = .true. elsewhere - pl%ldiscard(:) = .true. - pl%lmask(:) = .false. + pl%ldiscard(1:pl%nbody) = .true. + pl%lmask(1:pl%nbody) = .false. end where pl%nplm = count(pl%lmtiny(1:pl%nbody) .and. pl%lmask(1:pl%nbody))