From a7fb1c76168beb1827ebc1cc08dcc8b1c0b26e6f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 11 Dec 2022 17:21:35 -0500 Subject: [PATCH] I am once again restructuring. This time to try to get all of the storage objects to have a consistent structure and set of methods --- src/discard/discard.f90 | 4 +- src/encounter/encounter_io.f90 | 171 +++++++++++-------- src/encounter/encounter_util.f90 | 262 ++++++++++++++++++++++++++++++ src/fraggle/fraggle_io.f90 | 100 ++++++------ src/io/io.f90 | 101 ++++++------ src/main/swiftest_driver.f90 | 153 ++++++++--------- src/modules/encounter_classes.f90 | 55 +++++-- src/modules/fraggle_classes.f90 | 6 +- src/modules/swiftest_classes.f90 | 13 +- src/modules/symba_classes.f90 | 29 +--- src/netcdf/netcdf.f90 | 82 +++++----- src/setup/setup.f90 | 22 ++- src/symba/symba_collision.f90 | 6 +- src/symba/symba_io.f90 | 49 +----- src/symba/symba_step.f90 | 47 +++--- src/symba/symba_util.f90 | 247 +--------------------------- src/util/util_snapshot.f90 | 9 +- 17 files changed, 690 insertions(+), 666 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 1e0be68bd..72782df84 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 (ltp_discards) call tp_discards%write_info(param%nc, param) - if (lpl_discards) call pl_discards%write_info(param%nc, param) + if (ltp_discards) call tp_discards%write_info(param%system_history%nc, param) + if (lpl_discards) call pl_discards%write_info(param%system_history%nc, 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 index 0f821f49b..350f2f1d2 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -12,35 +12,49 @@ contains - module subroutine encounter_io_dump_collision_storage(self, param) + module subroutine encounter_io_dump_collision(self, param) !! author: David A. Minton !! !! Dumps the time history of an encounter to file. implicit none ! Arguments - class(collision_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(collision_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - select type(snapshot => self%frame(i)%item) - class is (fraggle_collision_snapshot) - param%ioutput = i - call snapshot%write_frame(self%nc,param) - end select - else - exit + select type(nc => self%nc) + class is (fraggle_io_parameters) + if (self%iframe > 0) then + nc%file_number = nc%file_number + 1 + nc%event_dimsize = self%iframe + call self%make_index_map() + write(nc%file_name, '("collision_",I0.6,".nc")') nc%file_number + call nc%initialize(param) + + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) then + select type(snapshot => self%frame(i)%item) + class is (fraggle_collision_snapshot) + param%ioutput = i + call snapshot%write_frame(nc,param) + end select + else + exit + end if + end do + + call nc%close() + call self%reset() end if - end do + end select return - end subroutine encounter_io_dump_collision_storage + end subroutine encounter_io_dump_collision - module subroutine encounter_io_dump_storage(self, param) - !! author: David A. Minton + module subroutine encounter_io_dump_encounter(self, param) + ! author: David A. Minton !! !! Dumps the time history of an encounter to file. implicit none @@ -50,20 +64,34 @@ module subroutine encounter_io_dump_storage(self, param) ! Internals integer(I4B) :: i - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - select type(snapshot => self%frame(i)%item) - class is (encounter_snapshot) - param%ioutput = self%tmap(i) - call snapshot%write_frame(self%nc,param) - end select - else - exit + select type(nc => self%nc) + class is (encounter_io_parameters) + if (self%iframe > 0) then + ! Create and save the output files for this encounter and fragmentation + nc%file_number = nc%file_number + 1 + call self%make_index_map() + write(nc%file_name, '("encounter_",I0.6,".nc")') nc%file_number + call nc%initialize(param) + + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) then + select type(snapshot => self%frame(i)%item) + class is (encounter_snapshot) + param%ioutput = self%tmap(i) + call snapshot%write_frame(nc,param) + end select + else + exit + end if + end do + + call nc%close() + call self%reset() end if - end do + end select return - end subroutine encounter_io_dump_storage + end subroutine encounter_io_dump_encounter module subroutine encounter_io_initialize(self, param) @@ -170,56 +198,59 @@ module subroutine encounter_io_write_frame(self, nc, param) use netcdf implicit none ! Arguments - class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp character(len=:), allocatable :: charstring tslot = param%ioutput associate(pl => self%pl, tp => self%tp) + select type (nc) + class is (encounter_io_parameters) - call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" ) - - call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) - call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame nf90_put_var pl loop_varid" ) - - npl = pl%nbody - do i = 1, npl - idslot = pl%id(i) + 1 - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl vh_varid" ) - call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid" ) - - if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid" ) - - if (param%lrotation) then - call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl Ip_varid" ) - call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rotx_varid" ) - end if - - charstring = trim(adjustl(pl%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var pl name_varid" ) - charstring = trim(adjustl(pl%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var pl particle_type_varid" ) - end do - - ntp = tp%nbody - do i = 1, ntp - idslot = tp%id(i) + 1 - call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var tp id_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp vh_varid" ) - - charstring = trim(adjustl(tp%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var tp name_varid" ) - charstring = trim(adjustl(tp%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var tp particle_type_varid" ) - end do - - call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" ) + + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame nf90_put_var pl loop_varid" ) + + npl = pl%nbody + do i = 1, npl + idslot = pl%id(i) + 1 + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid" ) + + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid" ) + + if (param%lrotation) then + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rotx_varid" ) + end if + + charstring = trim(adjustl(pl%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var pl name_varid" ) + charstring = trim(adjustl(pl%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var pl particle_type_varid" ) + end do + + ntp = tp%nbody + do i = 1, ntp + idslot = tp%id(i) + 1 + call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var tp id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp vh_varid" ) + + charstring = trim(adjustl(tp%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var tp name_varid" ) + charstring = trim(adjustl(tp%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var tp particle_type_varid" ) + end do + + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + end select end associate return diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index ef89b3cf3..45766a887 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -320,4 +320,266 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru return end subroutine encounter_util_spill_list + + + subroutine encounter_util_save_collision(param, snapshot) + !! author: David A. Minton + !! + !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. + implicit none + ! Arguments + class(symba_parameters), intent(inout) :: param !! SyMBA parameter object + class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object + ! Internals + type(collision_storage(nframes=:)), allocatable :: tmp + integer(I4B) :: i, nnew, nold, nbig + + ! Advance the snapshot frame counter + param%collision_history%iframe = param%collision_history%iframe + 1 + + ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 + nnew = param%collision_history%iframe + nold = param%collision_history%nframes + + if (nnew > nold) then + nbig = nold + do while (nbig < nnew) + nbig = nbig * 2 + end do + allocate(collision_storage(nbig) :: tmp) + tmp%iframe = param%collision_history%iframe + call move_alloc(param%collision_history%nc, tmp%nc) + + do i = 1, nold + if (allocated(param%collision_history%frame(i)%item)) call move_alloc(param%collision_history%frame(i)%item, tmp%frame(i)%item) + end do + deallocate(param%collision_history) + call move_alloc(tmp,param%collision_history) + nnew = nbig + end if + + param%collision_history%frame(nnew) = snapshot + + return + end subroutine encounter_util_save_collision + + + subroutine encounter_util_save_encounter(param, snapshot, t) + !! author: David A. Minton + !! + !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. + implicit none + ! Arguments + class(symba_parameters), intent(inout) :: param !! SyMBA parameter object + class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object + real(DP), intent(in) :: t !! The time of the snapshot + ! Internals + type(encounter_storage(nframes=:)), allocatable :: tmp + integer(I4B) :: i, nnew, nold, nbig + + ! Advance the snapshot frame counter + param%encounter_history%iframe = param%encounter_history%iframe + 1 + + ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 + nnew = param%encounter_history%iframe + nold = param%encounter_history%nframes + + if (nnew > nold) then + nbig = nold + do while (nbig < nnew) + nbig = nbig * 2 + end do + allocate(encounter_storage(nbig) :: tmp) + tmp%iframe = param%encounter_history%iframe + call move_alloc(param%encounter_history%nc, tmp%nc) + + do i = 1, nold + if (allocated(param%encounter_history%frame(i)%item)) call move_alloc(param%encounter_history%frame(i)%item, tmp%frame(i)%item) + end do + deallocate(param%encounter_history) + call move_alloc(tmp,param%encounter_history) + nnew = nbig + end if + + ! Find out which time slot this belongs in by searching for an existing slot + ! with the same value of time or the first available one + param%encounter_history%frame(nnew) = snapshot + + return + end subroutine encounter_util_save_encounter + + + module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) + !! author: David A. Minton + !! + !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories + !! can be played back through the encounter + implicit none + ! Internals + class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + ! Arguments + class(fraggle_collision_snapshot), allocatable :: snapshot + type(symba_pl) :: pl + character(len=:), allocatable :: stage + integer(I4B) :: i,j + + + if (present(arg)) then + stage = arg + else + stage = "" + end if + + select type (system) + class is (symba_nbody_system) + + select case(stage) + case("before") + ! Saves the states of the bodies involved in the collision before the collision is resolved + associate (idx => system%colliders%idx, ncoll => system%colliders%ncoll) + call pl%setup(ncoll, param) + pl%id(:) = system%pl%id(idx(:)) + pl%Gmass(:) = system%pl%Gmass(idx(:)) + pl%radius(:) = system%pl%radius(idx(:)) + pl%rot(:,:) = system%pl%rot(:,idx(:)) + pl%Ip(:,:) = system%pl%Ip(:,idx(:)) + pl%rh(:,:) = system%pl%rh(:,idx(:)) + pl%vh(:,:) = system%pl%vh(:,idx(:)) + pl%info(:) = system%pl%info(idx(:)) + !end select + allocate(system%colliders%pl, source=pl) + end associate + case("after") + allocate(fraggle_collision_snapshot :: snapshot) + allocate(snapshot%colliders, source=system%colliders) + allocate(snapshot%fragments, source=system%fragments) + select type (param) + class is (symba_parameters) + call encounter_util_save_collision(param,snapshot) + end select + case default + write(*,*) "encounter_util_snapshot_collision requies either 'before' or 'after' passed to 'arg'" + end select + + end select + + return + end subroutine encounter_util_snapshot_collision + + + module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) + !! author: David A. Minton + !! + !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories + !! can be played back through the encounter + implicit none + ! Internals + class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + ! Arguments + class(encounter_snapshot), allocatable :: snapshot + integer(I4B) :: i, npl_snap, ntp_snap + + if (.not.present(t)) then + write(*,*) "encounter_util_snapshot_encounter requires `t` to be passed" + end if + select type (system) + class is (symba_nbody_system) + select type(pl => system%pl) + class is (symba_pl) + select type (tp => system%tp) + class is (symba_tp) + associate(npl => pl%nbody, ntp => tp%nbody) + + allocate(encounter_snapshot :: snapshot) + snapshot%t = t + snapshot%iloop = param%iloop + + if (npl + ntp == 0) return + npl_snap = npl + ntp_snap = ntp + + allocate(snapshot%pl, mold=pl) + allocate(snapshot%tp, mold=tp) + select type(pl_snap => snapshot%pl) + class is (symba_pl) + if (npl > 0) then + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec + npl_snap = count(pl%lmask(1:npl)) + end if + if (ntp > 0) then + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec + ntp_snap = count(tp%lmask(1:ntp)) + end if + pl_snap%nbody = npl_snap + end select + + select type(pl_snap => snapshot%pl) + class is (symba_pl) + ! Take snapshot of the currently encountering massive bodies + if (npl_snap > 0) then + call pl_snap%setup(npl_snap, param) + pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) + pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) + pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) + pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) + do i = 1, NDIM + pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) + pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) + end do + if (param%lclose) then + pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) + end if + + if (param%lrotation) then + do i = 1, NDIM + pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) + pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) + end do + end if + call pl_snap%sort("id", ascending=.true.) + end if + end select + + select type(tp_snap => snapshot%tp) + class is (symba_tp) + ! Take snapshot of the currently encountering test particles + tp_snap%nbody = ntp_snap + if (ntp_snap > 0) then + call tp_snap%setup(ntp_snap, param) + tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) + tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) + do i = 1, NDIM + tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) + tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) + end do + end if + end select + end associate + ! Save the snapshot + select type(param) + class is (symba_parameters) + param%encounter_history%nid = param%encounter_history%nid + ntp_snap + npl_snap + call encounter_util_save_encounter(param,snapshot,t) + end select + end select + end select + end select + + return + end subroutine encounter_util_snapshot_encounter + end submodule s_encounter_util \ No newline at end of file diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 4f5fec19b..d4a8d3d9e 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -155,64 +155,66 @@ module subroutine fraggle_io_write_frame(self, nc, param) use netcdf implicit none ! Arguments - class(fraggle_collision_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(fraggle_collision_snapshot), intent(in) :: self !! Swiftest encounter structure + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, eslot, idslot, old_mode, npl, stage character(len=:), allocatable :: charstring class(swiftest_pl), allocatable :: pl associate(colliders => self%colliders, fragments => self%fragments) - eslot = param%ioutput select type(nc) - class is (fraggle_io_parameters) - call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "fraggle_io_write_frame nf90_set_fill" ) - - call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "fraggle_io_write_frame nf90_put_var time_varid" ) - call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "fraggle_io_write_frame nf90_put_varloop_varid" ) - - charstring = trim(adjustl(REGIME_NAMES(fragments%regime))) - call check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var regime_varid" ) - call check( nf90_put_var(nc%id, nc%Qloss_varid, fragments%Qloss, start=[eslot] ), "fraggle_io_write_frame nf90_put_var Qloss_varid" ) - - do stage = 1,2 - if (allocated(pl)) deallocate(pl) - select case(stage) - case(1) - allocate(pl, source=colliders%pl) - case(2) - allocate(pl, source=fragments%pl) - end select - npl = pl%nbody - do i = 1, npl - idslot = pl%id(i) + 1 - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "fraggle_io_write_frame nf90_put_var id_varid" ) - charstring = trim(adjustl(pl%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var name_varid" ) - charstring = trim(adjustl(pl%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], count=[len(charstring), 1, 1]), "fraggle_io_write_frame nf90_put_var particle_type_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var vh_varid" ) - call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "fraggle_io_write_frame nf90_put_var Gmass_varid" ) - call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "fraggle_io_write_frame nf90_put_var radius_varid" ) - call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var Ip_varid" ) - call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rotx_varid" ) + class is (encounter_io_parameters) + eslot = param%ioutput + select type(nc) + class is (fraggle_io_parameters) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "fraggle_io_write_frame nf90_set_fill" ) + + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "fraggle_io_write_frame nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "fraggle_io_write_frame nf90_put_varloop_varid" ) + + charstring = trim(adjustl(REGIME_NAMES(fragments%regime))) + call check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var regime_varid" ) + call check( nf90_put_var(nc%id, nc%Qloss_varid, fragments%Qloss, start=[eslot] ), "fraggle_io_write_frame nf90_put_var Qloss_varid" ) + + do stage = 1,2 + if (allocated(pl)) deallocate(pl) + select case(stage) + case(1) + allocate(pl, source=colliders%pl) + case(2) + allocate(pl, source=fragments%pl) + end select + npl = pl%nbody + do i = 1, npl + idslot = pl%id(i) + 1 + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "fraggle_io_write_frame nf90_put_var id_varid" ) + charstring = trim(adjustl(pl%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var name_varid" ) + charstring = trim(adjustl(pl%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], count=[len(charstring), 1, 1]), "fraggle_io_write_frame nf90_put_var particle_type_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "fraggle_io_write_frame nf90_put_var Gmass_varid" ) + call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "fraggle_io_write_frame nf90_put_var radius_varid" ) + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rotx_varid" ) + end do end do - end do - call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid before" ) - call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid after" ) - call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid before" ) - call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid after" ) - call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid before" ) - call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid after" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid before" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) - - call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid after" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid before" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid after" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid after" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid before" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) + + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + end select end select - - end associate + end associate return end subroutine fraggle_io_write_frame diff --git a/src/io/io.f90 b/src/io/io.f90 index 3d694e027..af4b30442 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -130,7 +130,7 @@ module subroutine io_conservation_report(self, param, lterminal) "; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, & "; DM/M0 = ", ES12.5)' - associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit) + associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit, nc => param%system_history%nc) call pl%vb2vh(cb) call pl%xh2xb(cb) @@ -177,8 +177,8 @@ 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 - call self%write_frame(param%nc, param) - call param%nc%close() + call self%write_frame(nc, param) + call nc%close() call util_exit(FAILURE) end if end if @@ -246,20 +246,22 @@ module subroutine io_dump_system(self, param) dump_param%out_stat = 'APPEND' dump_param%in_type = "NETCDF_DOUBLE" dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) - dump_param%nc%id_chunk = self%pl%nbody + self%tp%nbody - dump_param%nc%time_chunk = 1 - dump_param%tstart = self%t - - call dump_param%dump(param_file_name) - - dump_param%out_form = "XV" - dump_param%nc%file_name = trim(adjustl(DUMP_NC_FILE(idx))) - dump_param%ioutput = 1 - call dump_param%nc%initialize(dump_param) - call self%write_frame(dump_param%nc, dump_param) - call dump_param%nc%close() - ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) - call param%nc%flush(param) + associate(nc => dump_param%system_history%nc) + nc%id_chunk = self%pl%nbody + self%tp%nbody + nc%time_chunk = 1 + dump_param%tstart = self%t + + call dump_param%dump(param_file_name) + + dump_param%out_form = "XV" + nc%file_name = trim(adjustl(DUMP_NC_FILE(idx))) + dump_param%ioutput = 1 + call nc%initialize(dump_param) + call self%write_frame(nc, dump_param) + call nc%close() + ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) + call nc%flush(param) + end associate idx = idx + 1 if (idx > NDUMPFILES) idx = 1 @@ -267,14 +269,13 @@ module subroutine io_dump_system(self, param) ! Dump the encounter history if necessary select type(param) class is (symba_parameters) - if (param%lencounter_save) then - select type(self) - class is (symba_nbody_system) - call self%dump_encounter(param) - end select - end if + call param%encounter_history%dump(param) + call param%collision_history%dump(param) end select + ! Dump the system history to file + call param%system_history%dump(param) + return end subroutine io_dump_system @@ -1314,13 +1315,13 @@ module subroutine io_read_in_system(self, param) self%Euntracked = param%Euntracked else allocate(tmp_param, source=param) - tmp_param%nc%file_name = param%in_netcdf + tmp_param%system_history%nc%file_name = 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%nc, tmp_param) + ierr = self%read_frame(tmp_param%system_history%nc, tmp_param) deallocate(tmp_param) if (ierr /=0) call util_exit(FAILURE) end if @@ -1549,32 +1550,34 @@ module subroutine io_write_frame_system(self, param) character(len=STRMAX) :: errmsg logical :: fileExists - param%nc%id_chunk = self%pl%nbody + self%tp%nbody - param%nc%time_chunk = max(param%dump_cadence / param%istep_out, 1) - param%nc%file_name = param%outfile - if (lfirst) then - inquire(file=param%outfile, exist=fileExists) - - select case(param%out_stat) - case('APPEND') - if (.not.fileExists) then - errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" - goto 667 - end if - case('NEW') - if (fileExists) then - errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" - goto 667 - end if - call param%nc%initialize(param) - case('REPLACE', 'UNKNOWN') - call param%nc%initialize(param) - end select + associate (nc => param%system_history%nc, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) + nc%id_chunk = npl + ntp + nc%time_chunk = max(param%dump_cadence / param%istep_out, 1) + nc%file_name = param%outfile + if (lfirst) then + inquire(file=param%outfile, exist=fileExists) + + select case(param%out_stat) + case('APPEND') + if (.not.fileExists) then + errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" + goto 667 + end if + case('NEW') + if (fileExists) then + errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" + goto 667 + end if + call nc%initialize(param) + case('REPLACE', 'UNKNOWN') + call nc%initialize(param) + end select - lfirst = .false. - end if + lfirst = .false. + end if - call self%write_frame(param%nc, param) + call self%write_frame(nc, param) + end associate return diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 803ca6849..cab6d4aca 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -21,7 +21,7 @@ program swiftest_driver class(swiftest_nbody_system), allocatable :: system !! Polymorphic object containing the nbody system to be integrated class(swiftest_parameters), allocatable :: param !! Run configuration parameters character(len=:), allocatable :: integrator !! Integrator type code (see swiftest_globals for symbolic names) - character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters + character(len=:), allocatable :: param_file_name !! Name of the file containing user-defined parameters character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" integer(I8B) :: istart !! Starting index for loop counter integer(I8B) :: nloops !! Number of steps to take in the simulation @@ -38,7 +38,7 @@ program swiftest_driver character(len=64) :: pbarmessage character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' - type(swiftest_storage(nframes=:)), allocatable :: system_history + !type(swiftest_storage(nframes=:)), allocatable :: system_history call io_get_args(integrator, param_file_name, display_style) @@ -73,7 +73,6 @@ program swiftest_driver ! Set up system storage for intermittent file dumps if (dump_cadence == 0) dump_cadence = ceiling(nloops / (1.0_DP * istep_out), kind=I8B) - allocate(swiftest_storage(dump_cadence) :: system_history) ! Construct the main n-body system using the user-input integrator to choose the type of system call setup_construct_system(system, param) @@ -88,84 +87,86 @@ program swiftest_driver call system%initialize(param) - ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. - if (param%lrestart) then - if (param%lenergy) call system%conservation_report(param, lterminal=.true.) - else - if (param%lenergy) call system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum - call system_history%take_snapshot(system) - call system_history%dump(param) - end if - call system%dump(param) - - write(display_unit, *) " *************** Main Loop *************** " - - if (display_style == "PROGRESS") then - call pbar%reset(nloops) - write(pbarmessage,fmt=pbarfmt) t0, tstop - call pbar%update(1,message=pbarmessage) - else if (display_style == "COMPACT") then - write(*,*) "SWIFTEST START " // param%integrator - call system%compact_output(param,integration_timer) - end if - - iout = 0 - idump = 0 - system%t = tstart - do iloop = istart, nloops - !> Step the system forward in time - call integration_timer%start() - call system%step(param, system%t, dt) - call integration_timer%stop() - - system%t = t0 + iloop * dt - - !> Evaluate any discards or collisional outcomes - call system%discard(param) - if (display_style == "PROGRESS") call pbar%update(iloop) - - !> If the loop counter is at the output cadence value, append the data file with a single frame - if (istep_out > 0) then - iout = iout + 1 - if (iout == istep_out) then - iout = 0 - idump = idump + 1 - call system_history%take_snapshot(system) - - if (idump == dump_cadence) then - idump = 0 - call system%dump(param) - call system_history%dump(param) - end if + associate (system_history => param%system_history) + ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. + if (param%lrestart) then + if (param%lenergy) call system%conservation_report(param, lterminal=.true.) + else + if (param%lenergy) call system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum + call system_history%take_snapshot(param,system) + call system_history%dump(param) + end if + call system%dump(param) - tfrac = (system%t - t0) / (tstop - t0) - - select type(pl => system%pl) - class is (symba_pl) - write(display_unit, symbastatfmt) system%t, tfrac, pl%nplm, pl%nbody, system%tp%nbody - class default - write(display_unit, statusfmt) system%t, tfrac, pl%nbody, system%tp%nbody - end select - if (param%lenergy) call system%conservation_report(param, lterminal=.true.) - call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out) - - if (display_style == "PROGRESS") then - write(pbarmessage,fmt=pbarfmt) system%t, tstop - call pbar%update(1,message=pbarmessage) - else if (display_style == "COMPACT") then - call system%compact_output(param,integration_timer) - end if + write(display_unit, *) " *************** Main Loop *************** " - call integration_timer%reset() + if (display_style == "PROGRESS") then + call pbar%reset(nloops) + write(pbarmessage,fmt=pbarfmt) t0, tstop + call pbar%update(1,message=pbarmessage) + else if (display_style == "COMPACT") then + write(*,*) "SWIFTEST START " // param%integrator + call system%compact_output(param,integration_timer) + end if + iout = 0 + idump = 0 + system%t = tstart + do iloop = istart, nloops + !> Step the system forward in time + call integration_timer%start() + call system%step(param, system%t, dt) + call integration_timer%stop() + + system%t = t0 + iloop * dt + + !> Evaluate any discards or collisional outcomes + call system%discard(param) + if (display_style == "PROGRESS") call pbar%update(iloop) + + !> If the loop counter is at the output cadence value, append the data file with a single frame + if (istep_out > 0) then + iout = iout + 1 + if (iout == istep_out) then + iout = 0 + idump = idump + 1 + call system_history%take_snapshot(param,system) + + if (idump == dump_cadence) then + idump = 0 + call system%dump(param) + + end if + + tfrac = (system%t - t0) / (tstop - t0) + + select type(pl => system%pl) + class is (symba_pl) + write(display_unit, symbastatfmt) system%t, tfrac, pl%nplm, pl%nbody, system%tp%nbody + class default + write(display_unit, statusfmt) system%t, tfrac, pl%nbody, system%tp%nbody + end select + if (param%lenergy) call system%conservation_report(param, lterminal=.true.) + call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out) + + if (display_style == "PROGRESS") then + write(pbarmessage,fmt=pbarfmt) system%t, tstop + call pbar%update(1,message=pbarmessage) + else if (display_style == "COMPACT") then + call system%compact_output(param,integration_timer) + end if + + call integration_timer%reset() + + end if end if - end if - end do - ! Dump any remaining history if it exists - call system%dump(param) - call system_history%dump(param) - if (display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // param%integrator + end do + ! Dump any remaining history if it exists + call system%dump(param) + call system_history%dump(param) + if (display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // param%integrator + end associate end associate call util_exit(SUCCESS) diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index b1b7bb079..d339b6660 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -64,24 +64,24 @@ module encounter_classes procedure :: initialize => encounter_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object end type encounter_io_parameters - !> A class that that is used to store simulation history data between file output - type, extends(swiftest_storage) :: encounter_storage - class(encounter_io_parameters), allocatable :: nc !! NetCDF parameter object containing the details about the file attached to this storage object - contains - procedure :: dump => encounter_io_dump_storage !! Dumps contents of encounter history to file - procedure :: make_index_map => encounter_util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - final :: encounter_util_final_storage - end type encounter_storage - !> A class that that is used to store simulation history data between file output type, extends(swiftest_storage) :: collision_storage - class(encounter_io_parameters), allocatable :: nc !! NetCDF parameter object containing the details about the file attached to this storage object contains - procedure :: dump => encounter_io_dump_collision_storage !! Dumps contents of encounter history to file + procedure :: dump => encounter_io_dump_collision !! Dumps contents of encounter history to file procedure :: make_index_map => encounter_util_index_map_collision_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + procedure :: take_snapshot => encounter_util_snapshot_collision !! Take a minimal snapshot of the system through an encounter final :: encounter_util_final_collision_storage end type collision_storage + !> A class that that is used to store simulation history data between file output + type, extends(swiftest_storage) :: encounter_storage + contains + procedure :: dump => encounter_io_dump_encounter !! Dumps contents of encounter history to file + procedure :: make_index_map => encounter_util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + procedure :: take_snapshot => encounter_util_snapshot_encounter !! Take a minimal snapshot of the system through an encounter + final :: encounter_util_final_storage + end type encounter_storage + type encounter_bounding_box_1D integer(I4B) :: n !! Number of bodies with extents integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices (value > n indicates an ending index) @@ -214,17 +214,17 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_single_list - module subroutine encounter_io_dump_collision_storage(self, param) + module subroutine encounter_io_dump_collision(self, param) implicit none class(collision_storage(*)), intent(inout) :: self !! Collision storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_dump_collision_storage + end subroutine encounter_io_dump_collision - module subroutine encounter_io_dump_storage(self, param) + module subroutine encounter_io_dump_encounter(self, param) implicit none class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_dump_storage + end subroutine encounter_io_dump_encounter module subroutine encounter_io_initialize(self, param) implicit none @@ -234,9 +234,9 @@ end subroutine encounter_io_initialize module subroutine encounter_io_write_frame(self, nc, param) implicit none - class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine encounter_io_write_frame module subroutine encounter_setup_aabb(self, n, n_last) @@ -316,6 +316,25 @@ module subroutine encounter_util_resize_list(self, nnew) integer(I8B), intent(in) :: nnew !! New size of list needed end subroutine encounter_util_resize_list + module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) + implicit none + class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + end subroutine encounter_util_snapshot_collision + + module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) + implicit none + class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + end subroutine encounter_util_snapshot_encounter + + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 8d8e1c33e..2bd120b83 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -160,9 +160,9 @@ end subroutine fraggle_io_initialize_output module subroutine fraggle_io_write_frame(self, nc, param) implicit none - class(fraggle_collision_snapshot), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(fraggle_collision_snapshot), intent(in) :: self !! Swiftest encounter structure + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine fraggle_io_write_frame module subroutine fraggle_io_log_pl(pl, param) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4894e305d..b89e9eee8 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -162,6 +162,7 @@ module swiftest_classes integer(I4B) :: nt !! Number of unique time values in all saved snapshots real(DP), dimension(:), allocatable :: tvals !! The set of unique time values contained in the snapshots integer(I4B), dimension(:), allocatable :: tmap !! The t value -> index map + class(netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: dump => io_dump_storage !! Dumps storage object contents to file procedure :: make_index_map => util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id @@ -255,8 +256,7 @@ module swiftest_classes logical :: lgr = .false. !! Turn on GR logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect logical :: lyorp = .false. !! Turn on YORP effect - - type(netcdf_parameters) :: nc !! Object containing NetCDF parameters + type(swiftest_storage(nframes=:)), allocatable :: system_history contains procedure :: reader => io_param_reader procedure :: writer => io_param_writer @@ -1692,10 +1692,13 @@ module subroutine util_set_rhill_approximate(self,cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine util_set_rhill_approximate - module subroutine util_snapshot_system(self, system) + module subroutine util_snapshot_system(self, param, system, t, arg) implicit none - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object - class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object to store + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) end subroutine util_snapshot_system end interface diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index ff2597b5f..3dbcfb90f 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -33,6 +33,8 @@ module symba_classes character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved character(STRMAX) :: collision_save = "NONE" !! Indicate if and how fragmentation data should be saved logical :: lencounter_save = .false. !! Turns on encounter saving + type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file + type(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file contains procedure :: reader => symba_io_param_reader procedure :: writer => symba_io_param_writer @@ -191,8 +193,6 @@ module symba_classes integer(I4B) :: irec !! System recursion level class(fraggle_colliders), allocatable :: colliders !! Fraggle colliders object class(fraggle_fragments), allocatable :: fragments !! Fraggle fragmentation system object - type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file - type(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps @@ -201,9 +201,6 @@ module symba_classes procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current system level procedure :: recursive_step => symba_step_recur_system !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary procedure :: reset => symba_step_reset_system !! Resets pl, tp,and encounter structures at the start of a new step - procedure :: encounter_snap => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter - procedure :: collision_snap => symba_util_take_collision_snapshot !! Take a minimal snapshot of the system before and after a collision - procedure :: dump_encounter => symba_io_dump_encounter !! Saves the encounter and/or fragmentation data to file(s) final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system @@ -373,22 +370,6 @@ module subroutine symba_util_set_renc(self, scale) integer(I4B), intent(in) :: scale !! Current recursion depth end subroutine symba_util_set_renc - module subroutine symba_util_take_collision_snapshot(self, param, t, stage) - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - character(*), intent(in) :: stage !! Either before or afte - end subroutine symba_util_take_collision_snapshot - - module subroutine symba_util_take_encounter_snapshot(self, param, t) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - end subroutine symba_util_take_encounter_snapshot - module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none class(symba_parameters), intent(inout) :: self !! Current run configuration parameters with SyMBA additionss @@ -411,12 +392,6 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine symba_io_param_writer - module subroutine symba_io_dump_encounter(self, param) - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_io_dump_encounter - module subroutine symba_io_write_discard(self, param) use swiftest_classes, only : swiftest_parameters implicit none diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 9bc945227..42e2a2ea6 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -80,63 +80,65 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) real(DP), dimension(NDIM) :: rot0, Ip0, Lnow real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig - call param%nc%open(param) - call check( nf90_inquire_dimension(param%nc%id, param%nc%time_dimid, len=itmax), "netcdf_get_old_t_final_system time_dimid" ) - call check( nf90_inquire_dimension(param%nc%id, param%nc%id_dimid, len=idmax), "netcdf_get_old_t_final_system id_dimid" ) - allocate(vals(idmax)) - call check( nf90_get_var(param%nc%id, param%nc%time_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system time_varid" ) + associate (nc => param%system_history%nc) + call nc%open(param) + call check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=itmax), "netcdf_get_old_t_final_system time_dimid" ) + call check( nf90_inquire_dimension(nc%id, nc%id_dimid, len=idmax), "netcdf_get_old_t_final_system id_dimid" ) + allocate(vals(idmax)) + call check( nf90_get_var(nc%id, nc%time_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system time_varid" ) - !old_t_final = rtemp(1) - old_t_final = param%t0 ! For NetCDF it is safe to overwrite the final t value on a restart + !old_t_final = rtemp(1) + old_t_final = param%t0 ! For NetCDF it is safe to overwrite the final t value on a restart - if (param%lenergy) then - call check( nf90_get_var(param%nc%id, param%nc%KE_orb_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_orb_varid" ) - KE_orb_orig = rtemp(1) + if (param%lenergy) then + call check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_orb_varid" ) + KE_orb_orig = rtemp(1) - call check( nf90_get_var(param%nc%id, param%nc%KE_spin_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_spin_varid" ) - KE_spin_orig = rtemp(1) + call check( nf90_get_var(nc%id, nc%KE_spin_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system KE_spin_varid" ) + KE_spin_orig = rtemp(1) - call check( nf90_get_var(param%nc%id, param%nc%PE_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system PE_varid" ) - PE_orig = rtemp(1) + call check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[1], count=[1]), "netcdf_get_old_t_final_system PE_varid" ) + PE_orig = rtemp(1) - call check( nf90_get_var(param%nc%id, param%nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_get_old_t_final_system Ecollisions_varid" ) - call check( nf90_get_var(param%nc%id, param%nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_get_old_t_final_system Euntracked_varid" ) + call check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_get_old_t_final_system Ecollisions_varid" ) + call check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_get_old_t_final_system Euntracked_varid" ) - self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked + self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked - call check( nf90_get_var(param%nc%id, param%nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_orb_varid" ) - call check( nf90_get_var(param%nc%id, param%nc%L_spin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_spin_varid" ) - call check( nf90_get_var(param%nc%id, param%nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_escape_varid" ) + call check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_orb_varid" ) + call check( nf90_get_var(nc%id, nc%L_spin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_spin_varid" ) + call check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_escape_varid" ) - self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) + self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) - call check( nf90_get_var(param%nc%id, param%nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_get_old_t_final_system Gmass_varid" ) - call check( nf90_get_var(param%nc%id, param%nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_get_old_t_final_system GMescape_varid" ) - self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + self%GMescape + call check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_get_old_t_final_system Gmass_varid" ) + call check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_get_old_t_final_system GMescape_varid" ) + self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + self%GMescape - select type(cb => self%cb) - class is (symba_cb) - cb%GM0 = vals(1) - cb%dGM = cb%Gmass - cb%GM0 + select type(cb => self%cb) + class is (symba_cb) + cb%GM0 = vals(1) + cb%dGM = cb%Gmass - cb%GM0 - call check( nf90_get_var(param%nc%id, param%nc%radius_varid, rtemp, start=[1,1], count=[1,1]), "netcdf_get_old_t_final_system radius_varid" ) - cb%R0 = rtemp(1) + call check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,1], count=[1,1]), "netcdf_get_old_t_final_system radius_varid" ) + cb%R0 = rtemp(1) - if (param%lrotation) then + if (param%lrotation) then - call check( nf90_get_var(param%nc%id, param%nc%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system rot_varid" ) - call check( nf90_get_var(param%nc%id, param%nc%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system Ip_varid" ) + call check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system rot_varid" ) + call check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_get_old_t_final_system Ip_varid" ) - cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:) + cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:) - Lnow(:) = cb%Ip(3) * cb%Gmass * cb%radius**2 * cb%rot(:) - cb%dL(:) = Lnow(:) - cb%L0(:) - end if - end select + Lnow(:) = cb%Ip(3) * cb%Gmass * cb%radius**2 * cb%rot(:) + cb%dL(:) = Lnow(:) - cb%L0(:) + end if + end select - end if + end if - deallocate(vals) + deallocate(vals) + end associate return end function netcdf_get_old_t_final_system diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index a107179a1..64c63b390 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -21,6 +21,10 @@ module subroutine setup_construct_system(system, param) class(swiftest_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Swiftest parameters + allocate(swiftest_storage(param%dump_cadence) :: param%system_history) + allocate(netcdf_parameters :: param%system_history%nc) + call param%system_history%reset() + select case(param%integrator) case (BS) write(*,*) 'Bulirsch-Stoer integrator not yet enabled' @@ -72,8 +76,8 @@ module subroutine setup_construct_system(system, param) select type(param) class is (symba_parameters) if (param%lencounter_save) then - allocate(encounter_storage :: system%encounter_history) - associate (encounter_history => system%encounter_history) + allocate(encounter_storage :: param%encounter_history) + associate (encounter_history => param%encounter_history) allocate(encounter_io_parameters :: encounter_history%nc) call encounter_history%reset() select type(nc => encounter_history%nc) @@ -81,9 +85,9 @@ module subroutine setup_construct_system(system, param) nc%file_number = param%iloop / param%dump_cadence end select end associate - - allocate(collision_storage :: system%collision_history) - associate (collision_history => system%collision_history) + + allocate(collision_storage :: param%collision_history) + associate (collision_history => param%collision_history) allocate(fraggle_io_parameters :: collision_history%nc) call collision_history%reset() select type(nc => collision_history%nc) @@ -93,6 +97,8 @@ module subroutine setup_construct_system(system, param) end associate end if end select + + end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' @@ -101,6 +107,10 @@ module subroutine setup_construct_system(system, param) call util_exit(FAILURE) end select + + + + return end subroutine setup_construct_system @@ -116,7 +126,7 @@ module subroutine setup_finalize_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters associate(system => self) - call param%nc%close() + call param%system_history%nc%close() end associate return diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 95d0dcf4f..ab2962054 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -885,7 +885,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) logical :: lgoodcollision integer(I4B) :: i - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, t => system%t) + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, t => system%t, collision_history => param%collision_history) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) @@ -900,7 +900,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) call system%colliders%regime(system%fragments, system, param) - if (param%lencounter_save) call system%collision_snap(param, t, "before") + if (param%lencounter_save) call collision_history%take_snapshot(param,system, t, "before") select case (system%fragments%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) plplcollision_list%status(i) = symba_collision_casedisruption(system, param) @@ -912,7 +912,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) write(*,*) "Error in symba_collision, unrecognized collision regime" call util_exit(FAILURE) end select - if (param%lencounter_save) call system%collision_snap(param, t, "after") + if (param%lencounter_save) call collision_history%take_snapshot(param,system, t, "after") deallocate(system%colliders,system%fragments) end do end select diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index d1f36abca..18b56767e 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -12,51 +12,6 @@ contains - module subroutine symba_io_dump_encounter(self, param) - !! author: David A. Minton - !! - !! Saves the encounter and/or fragmentation data to file with the name of the current output interval number attached - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters - - - associate(encounter_history => self%encounter_history, num_enc_frames => self%encounter_history%iframe,& - collision_history => self%collision_history, num_coll_frames => self%collision_history%iframe) - - select type(nce => self%encounter_history%nc) - class is (encounter_io_parameters) - if (num_enc_frames > 0) then - ! Create and save the output files for this encounter and fragmentation - nce%file_number = nce%file_number + 1 - call encounter_history%make_index_map() - write(nce%file_name, '("encounter_",I0.6,".nc")') nce%file_number - call nce%initialize(param) - call encounter_history%dump(param) - call nce%close() - call encounter_history%reset() - end if - end select - - select type(ncc => self%collision_history%nc) - class is (fraggle_io_parameters) - if (num_coll_frames > 0) then - ncc%file_number = ncc%file_number + 1 - ncc%event_dimsize = num_coll_frames - call collision_history%make_index_map() - write(ncc%file_name, '("collision_",I0.6,".nc")') ncc%file_number - call ncc%initialize(param) - call collision_history%dump(param) - call ncc%close() - call collision_history%reset() - end if - end select - end associate - - return - end subroutine symba_io_dump_encounter - module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott @@ -246,12 +201,12 @@ module subroutine symba_io_write_discard(self, param) associate(pl => self%pl, npl => self%pl%nbody, pl_adds => self%pl_adds) - if (self%tp_discards%nbody > 0) call self%tp_discards%write_info(param%nc, param) + if (self%tp_discards%nbody > 0) call self%tp_discards%write_info(param%system_history%nc, param) select type(pl_discards => self%pl_discards) class is (symba_merger) if (pl_discards%nbody == 0) return - call pl_discards%write_info(param%nc, param) + call pl_discards%write_info(param%system_history%nc, param) end select end associate diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index dc303b4f7..54e2464d1 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -34,17 +34,19 @@ module subroutine symba_step_system(self, param, t, dt) class is (symba_tp) select type(param) class is (symba_parameters) - call self%reset(param) - lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) - if (lencounter) then - if (param%lencounter_save) call self%encounter_snap(param, t) - call self%interp(param, t, dt) - if (param%lencounter_save) call self%encounter_snap(param, t+dt) - else - self%irec = -1 - call helio_step_system(self, param, t, dt) - end if - param%lfirstkick = pl%lfirst + associate(encounter_history => param%encounter_history) + call self%reset(param) + lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) + if (lencounter) then + if (param%lencounter_save) call encounter_history%take_snapshot(param, self, t) + call self%interp(param, t, dt) + if (param%lencounter_save) call encounter_history%take_snapshot(param, self, t+dt) + else + self%irec = -1 + call helio_step_system(self, param, t, dt) + end if + param%lfirstkick = pl%lfirst + end associate end select end select end select @@ -180,14 +182,15 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) real(DP) :: dtl, dth logical :: lencounter - associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, & - lplpl_collision => self%plplenc_list%lcollision, lpltp_collision => self%pltpenc_list%lcollision) - select type(param) - class is (symba_parameters) - select type(pl => self%pl) - class is (symba_pl) - select type(tp => self%tp) - class is (symba_tp) + select type(param) + class is (symba_parameters) + select type(pl => self%pl) + class is (symba_pl) + select type(tp => self%tp) + class is (symba_tp) + associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, & + lplpl_collision => self%plplenc_list%lcollision, lpltp_collision => self%pltpenc_list%lcollision, & + encounter_history => param%encounter_history) system%irec = ireci dtl = param%dt / (NTENC**ireci) dth = 0.5_DP * dtl @@ -244,15 +247,15 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) end if - if (param%lencounter_save) call system%encounter_snap(param, t+dtl) + if (param%lencounter_save) call encounter_history%take_snapshot(param, self, t+dtl) call self%set_recur_levels(ireci) end do - end select + end associate end select end select - end associate + end select return end subroutine symba_step_recur_system diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 4381ce93b..9443d658c 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -667,7 +667,7 @@ module subroutine symba_util_rearray_pl(self, system, param) end where end select - call pl%write_info(param%nc, param) + call pl%write_info(param%system_history%nc, param) deallocate(ldump_mask) ! Reindex the new list of bodies @@ -868,100 +868,6 @@ module subroutine symba_util_resize_pl(self, nnew) return end subroutine symba_util_resize_pl - - subroutine symba_util_save_collision(system, snapshot) - !! author: David A. Minton - !! - !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an - !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment - !! Memory usage grows by a factor of 2 each time it fills up, but no more. - implicit none - ! Arguments - type(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object - ! Internals - type(collision_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nnew, nold, nbig - - ! Advance the snapshot frame counter - system%collision_history%iframe = system%collision_history%iframe + 1 - - ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = system%collision_history%iframe - nold = system%collision_history%nframes - - if (nnew > nold) then - nbig = nold - do while (nbig < nnew) - nbig = nbig * 2 - end do - allocate(collision_storage(nbig) :: tmp) - tmp%iframe = system%collision_history%iframe - call move_alloc(system%collision_history%nc, tmp%nc) - - do i = 1, nold - if (allocated(system%collision_history%frame(i)%item)) call move_alloc(system%collision_history%frame(i)%item, tmp%frame(i)%item) - end do - deallocate(system%collision_history) - call move_alloc(tmp,system%collision_history) - nnew = nbig - end if - - system%collision_history%frame(nnew) = snapshot - - return - end subroutine symba_util_save_collision - - - subroutine symba_util_save_encounter(system, snapshot, t) - !! author: David A. Minton - !! - !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an - !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment - !! Memory usage grows by a factor of 2 each time it fills up, but no more. - implicit none - ! Arguments - type(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object - real(DP), intent(in) :: t !! The time of the snapshot - ! Internals - type(encounter_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nnew, nold, nbig - - ! Advance the snapshot frame counter - system%encounter_history%iframe = system%encounter_history%iframe + 1 - - ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = system%encounter_history%iframe - nold = system%encounter_history%nframes - - if (nnew > nold) then - nbig = nold - do while (nbig < nnew) - nbig = nbig * 2 - end do - allocate(encounter_storage(nbig) :: tmp) - tmp%iframe = system%encounter_history%iframe - call move_alloc(system%encounter_history%nc, tmp%nc) - - do i = 1, nold - if (allocated(system%encounter_history%frame(i)%item)) call move_alloc(system%encounter_history%frame(i)%item, tmp%frame(i)%item) - end do - deallocate(system%encounter_history) - call move_alloc(tmp,system%encounter_history) - nnew = nbig - end if - - ! Find out which time slot this belongs in by searching for an existing slot - ! with the same value of time or the first available one - system%encounter_history%frame(nnew) = snapshot - - return - end subroutine symba_util_save_encounter - - module subroutine symba_util_resize_tp(self, nnew) !! author: David A. Minton !! @@ -1320,156 +1226,5 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) end subroutine symba_util_spill_tp - module subroutine symba_util_take_collision_snapshot(self, param, t, stage) - !! author: David A. Minton - !! - !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories - !! can be played back through the encounter - implicit none - ! Internals - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - character(*), intent(in) :: stage !! Either before or after - ! Arguments - class(fraggle_collision_snapshot), allocatable :: snapshot - type(symba_pl) :: pl - integer(I4B) :: i,j - - select case(stage) - case("before") - ! Saves the states of the bodies involved in the collision before the collision is resolved - associate (idx => self%colliders%idx, ncoll => self%colliders%ncoll) - !allocate(symba_pl :: self%colliders%pl) - !select type(pl => self%colliders%pl) - !class is (symba_pl) - call pl%setup(ncoll, param) - pl%id(:) = self%pl%id(idx(:)) - pl%Gmass(:) = self%pl%Gmass(idx(:)) - pl%radius(:) = self%pl%radius(idx(:)) - pl%rot(:,:) = self%pl%rot(:,idx(:)) - pl%Ip(:,:) = self%pl%Ip(:,idx(:)) - pl%rh(:,:) = self%pl%rh(:,idx(:)) - pl%vh(:,:) = self%pl%vh(:,idx(:)) - pl%info(:) = self%pl%info(idx(:)) - !end select - allocate(self%colliders%pl, source=pl) - end associate - case("after") - allocate(fraggle_collision_snapshot :: snapshot) - allocate(snapshot%colliders, source=self%colliders) - allocate(snapshot%fragments, source=self%fragments) - call symba_util_save_collision(self,snapshot) - end select - - return - end subroutine symba_util_take_collision_snapshot - - - module subroutine symba_util_take_encounter_snapshot(self, param, t) - !! author: David A. Minton - !! - !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories - !! can be played back through the encounter - implicit none - ! Internals - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - ! Arguments - class(encounter_snapshot), allocatable :: snapshot - integer(I4B) :: i, npl_snap, ntp_snap - - associate(npl => self%pl%nbody, ntp => self%tp%nbody) - - allocate(encounter_snapshot :: snapshot) - snapshot%t = t - snapshot%iloop = param%iloop - - if (npl + ntp == 0) return - npl_snap = npl - ntp_snap = ntp - - select type (pl => self%pl) - class is (symba_pl) - select type (tp => self%tp) - class is (symba_tp) - allocate(symba_pl :: snapshot%pl) - allocate(symba_tp :: snapshot%tp) - - select type(pl_snap => snapshot%pl) - class is (symba_pl) - select type(tp_snap => snapshot%tp) - class is (symba_tp) - - if (npl > 0) then - pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == self%irec - npl_snap = count(pl%lmask(1:npl)) - end if - if (ntp > 0) then - tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == self%irec - ntp_snap = count(tp%lmask(1:ntp)) - end if - pl_snap%nbody = npl_snap - - ! Take snapshot of the currently encountering massive bodies - if (npl_snap > 0) then - allocate(pl_snap%id(npl_snap)) - allocate(pl_snap%info(npl_snap)) - allocate(pl_snap%Gmass(npl_snap)) - - allocate(pl_snap%levelg(npl_snap)) - pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) - pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) - pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) - pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) - allocate(pl_snap%rh(NDIM,npl_snap)) - allocate(pl_snap%vh(NDIM,npl_snap)) - do i = 1, NDIM - pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) - pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) - end do - if (param%lclose) then - allocate(pl_snap%radius(npl_snap)) - pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) - end if - - if (param%lrotation) then - allocate(pl_snap%Ip(NDIM,npl_snap)) - allocate(pl_snap%rot(NDIM,npl_snap)) - do i = 1, NDIM - pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) - pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) - end do - end if - call pl_snap%sort("id", ascending=.true.) - end if - - ! Take snapshot of the currently encountering test particles - tp_snap%nbody = ntp_snap - if (ntp_snap > 0) then - allocate(tp_snap%id(ntp_snap)) - allocate(tp_snap%info(ntp_snap)) - tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) - tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) - allocate(tp_snap%rh(NDIM,ntp_snap)) - allocate(tp_snap%vh(NDIM,ntp_snap)) - do i = 1, NDIM - tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) - tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) - end do - end if - end select - end select - - ! Save the snapshot - self%encounter_history%nid = self%encounter_history%nid + ntp_snap + npl_snap - call symba_util_save_encounter(self,snapshot,t) - end select - end select - end associate - - return - end subroutine symba_util_take_encounter_snapshot end submodule s_symba_util diff --git a/src/util/util_snapshot.f90 b/src/util/util_snapshot.f90 index 1c67bc7f8..c3a98855b 100644 --- a/src/util/util_snapshot.f90 +++ b/src/util/util_snapshot.f90 @@ -11,14 +11,17 @@ use swiftest contains - module subroutine util_snapshot_system(self, system) + module subroutine util_snapshot_system(self, param, system, t, arg) !! author: David A. Minton !! !! Takes a snapshot of the system for later file storage implicit none ! Arguments - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object - class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object to store + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) self%iframe = self%iframe + 1 self%nt = self%iframe