Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Added energy file to parameter read. Started particle info io methods…
Browse files Browse the repository at this point in the history 
… in SyMBA
  • Loading branch information
@daminton
daminton committed Aug 10, 2021
1 parent 44df3c5 commit 88ba918
Show file tree
Hide file tree
Showing 6 changed files with 139 additions and 81 deletions.
18 changes: 11 additions & 7 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,11 @@ module subroutine io_conservation_report(self, param, lterminal)
Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig, &
Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), &
lfirst => param%lfirstenergy)
if (lfirst) then
if (param%out_stat == "OLD") then
open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "old", action = "write", position = "append")
else
open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "replace", action = "write")
if (param%energy_out /= "") then
if (lfirst .and. (param%out_stat /= "OLD")) then
open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write")
else
open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append")
write(EGYIU,EGYHEADER)
end if
end if
Expand All @@ -59,8 +59,10 @@ module subroutine io_conservation_report(self, param, lterminal)
lfirst = .false.
end if

write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
flush(EGYIU)
if (param%energy_out /= "") then
write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
close(EGYIU)
end if
if (.not.lfirst .and. lterminal) then
Lmag_now = norm2(Ltot_now)
Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig
Expand Down Expand Up @@ -422,6 +424,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
param%enc_out = param_value
case ("DISCARD_OUT")
param%discard_out = param_value
case ("ENERGY_OUT")
param%energy_out = param_value
case ("EXTRA_FORCE")
call io_toupper(param_value)
if (param_value == "YES" .or. param_value == 'T') param%lextra_force = .true.
Expand Down
2 changes: 1 addition & 1 deletion src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ module swiftest_classes
real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters
real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units
real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units
character(STRMAX) :: ennergy_out = "" !! Name of output energy and momentum report file
character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file

! Logical flags to turn on or off various features of the code
logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually)
Expand Down
77 changes: 33 additions & 44 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,33 +27,17 @@ module symba_classes
procedure :: writer => symba_io_param_writer
end type symba_parameters

!********************************************************************************************************************************
! symba_cb class definitions and method interfaces
!*******************************************************************************************************************************
!> SyMBA central body particle class
type, extends(helio_cb) :: symba_cb
real(DP) :: M0 = 0.0_DP !! Initial mass of the central body
real(DP) :: dM = 0.0_DP !! Change in mass of the central body
real(DP) :: R0 = 0.0_DP !! Initial radius of the central body
real(DP) :: dR = 0.0_DP !! Change in the radius of the central body
contains
end type symba_cb

!********************************************************************************************************************************
! symba_particle_info class definitions and method interfaces
!*******************************************************************************************************************************
!> Class definition for the particle origin information object. This object is used to track time, location, and collisional regime
!> of fragments produced in collisional events.
type, extends(swiftest_base) :: symba_particle_info
type :: symba_particle_info
sequence
character(len=32) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.)
real(DP) :: origin_time !! The time of the particle's formation
real(DP), dimension(NDIM) :: origin_xh !! The heliocentric distance vector at the time of the particle's formation
real(DP), dimension(NDIM) :: origin_vh !! The heliocentric velocity vector at the time of the particle's formation
contains
procedure :: dump => symba_io_dump_particle_info !! I/O routine for dumping particle info to file
procedure :: initialize => symba_io_initialize_particle_info !! I/O routine for reading in particle info data
procedure :: read_frame => symba_io_read_frame_info !! I/O routine for reading in a single frame of particle info
procedure :: write_frame => symba_io_write_frame_info !! I/O routine for writing out a single frame of particle info
end type symba_particle_info

!********************************************************************************************************************************
Expand All @@ -66,6 +50,19 @@ module symba_classes
integer(I4B), dimension(:), allocatable :: child !! Index of children particles
end type symba_kinship

!********************************************************************************************************************************
! symba_cb class definitions and method interfaces
!*******************************************************************************************************************************
!> SyMBA central body particle class
type, extends(helio_cb) :: symba_cb
real(DP) :: M0 = 0.0_DP !! Initial mass of the central body
real(DP) :: dM = 0.0_DP !! Change in mass of the central body
real(DP) :: R0 = 0.0_DP !! Initial radius of the central body
real(DP) :: dR = 0.0_DP !! Change in the radius of the central body
type(symba_particle_info) :: info
contains
end type symba_cb

!********************************************************************************************************************************
! symba_pl class definitions and method interfaces
!*******************************************************************************************************************************
Expand Down Expand Up @@ -118,6 +115,7 @@ module symba_classes
integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with planets this time step
integer(I4B), dimension(:), allocatable :: levelg !! level at which this particle should be moved
integer(I4B), dimension(:), allocatable :: levelm !! deepest encounter level achieved this time step
type(symba_particle_info), dimension(:), allocatable :: info
contains
procedure :: drift => symba_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body
Expand Down Expand Up @@ -327,12 +325,12 @@ module subroutine symba_io_write_discard(self, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine symba_io_write_discard

module subroutine symba_io_dump_particle_info(self, param, msg)
use swiftest_classes, only : swiftest_parameters
module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx)
implicit none
class(symba_particle_info), intent(inout) :: self !! Swiftest base object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), optional, intent(in) :: msg !! Message to display with dump operation
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions
integer(I4B), dimension(:), optional, intent(in) :: tpidx !! Array of test particle indices to append to the particle file
integer(I4B), dimension(:), optional, intent(in) :: plidx !! Array of massive body indices to append to the particle file
end subroutine symba_io_dump_particle_info

module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
Expand All @@ -357,22 +355,21 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms
character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0
end subroutine symba_io_param_writer

module subroutine symba_io_initialize_particle_info(self, param)
use swiftest_classes, only : swiftest_parameters
module subroutine symba_io_initialize_particle_info(system, param)
implicit none
class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions
end subroutine symba_io_initialize_particle_info

module subroutine symba_io_read_frame_info(self, iu, param, form, ierr)
use swiftest_classes, only : swiftest_parameters
implicit none
class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), intent(in) :: form !! Input format code ("XV" or "EL")
integer(I4B), intent(out) :: ierr !! Error code
end subroutine symba_io_read_frame_info
!module subroutine symba_io_read_frame_info(self, iu, param, form, ierr)
! use swiftest_classes, only : swiftest_parameters
! implicit none
! class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object
! integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
! class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! character(*), intent(in) :: form !! Input format code ("XV" or "EL")
! integer(I4B), intent(out) :: ierr !! Error code
!end subroutine symba_io_read_frame_info

module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
Expand Down Expand Up @@ -403,14 +400,6 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_pltpenc

module subroutine symba_io_write_frame_info(self, iu, param)
use swiftest_classes, only : swiftest_parameters
implicit none
class(symba_particle_info), intent(in) :: self !! SyMBA particle info object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine symba_io_write_frame_info

module subroutine symba_setup_initialize_system(self, param)
use swiftest_classes, only : swiftest_parameters
implicit none
Expand Down
116 changes: 87 additions & 29 deletions src/symba/symba_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,25 +2,91 @@
use swiftest
contains

module subroutine symba_io_dump_particle_info(self, param, msg)
module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx)
!! author: David A. Minton
!!
!! Dumps the particle information data to a file
!! Dumps the particle information data to a file.
!! Pass a list of array indices for test particles (tpidx) and/or massive bodies (plidx) to append
implicit none
class(symba_particle_info), intent(inout) :: self !! Swiftest base object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), optional, intent(in) :: msg !! Message to display with dump operation
! Arguments
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions
integer(I4B), dimension(:), optional, intent(in) :: tpidx !! Array of test particle indices to append to the particle file
integer(I4B), dimension(:), optional, intent(in) :: plidx !! Array of massive body indices to append to the particle file
! Internals
logical, save :: lfirst = .true.
integer(I4B), parameter :: iu = 22
integer(I4B) :: i, ierr

if (.not.present(tpidx) .and. .not.present(plidx)) return
if (lfirst) then
select case(param%out_stat)
case('APPEND')
open(unit = iu, file = param%particle_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr)
case('NEW', 'UNKNOWN', 'REPLACE')
open(unit = iu, file = param%particle_file, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr)
case default
write(*,*) 'Invalid status code',trim(adjustl(param%out_stat))
call util_exit(FAILURE)
end select
if (ierr /= 0) then
write(*, *) "Swiftest error:"
write(*, *) " particle output file already exists or cannot be accessed"
call util_exit(FAILURE)
end if

lfirst = .false.
else
open(unit = iu, file = param%particle_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr)
if (ierr /= 0) then
write(*, *) "Swiftest error:"
write(*, *) " unable to open binary output file for APPEND"
call util_exit(FAILURE)
end if
end if

if (present(plidx) .and. (system%pl%nbody > 0) .and. size(plidx) > 0) then
select type(pl => system%pl)
class is (symba_pl)
do i = 1, size(plidx)
write(iu) pl%id(plidx(i))
write(iu) pl%info(plidx(i))
end do
end select
end if

if (present(tpidx) .and. (system%tp%nbody > 0) .and. size(tpidx) > 0) then
select type(tp => system%tp)
class is (symba_tp)
do i = 1, size(tpidx)
write(iu) tp%id(tpidx(i))
write(iu) tp%info(tpidx(i))
end do
end select
end if

close(unit = iu, iostat = ierr)
if (ierr /= 0) then
write(*, *) "Swiftest error:"
write(*, *) " unable to close particle output file"
call util_exit(FAILURE)
end if

return
end subroutine symba_io_dump_particle_info


module subroutine symba_io_initialize_particle_info(self, param)
module subroutine symba_io_initialize_particle_info(system, param)
!! author: David A. Minton
!!
!! Initializes a particle info data structure, either starting a new one or reading one in
!! from a file if it is a restarted run
implicit none
class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Argumets
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions

return
end subroutine symba_io_initialize_particle_info


Expand Down Expand Up @@ -194,19 +260,19 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms
end subroutine symba_io_param_writer


module subroutine symba_io_read_frame_info(self, iu, param, form, ierr)
!! author: David A. Minton
!!
!! Reads a single frame of a particle info data from a file.
implicit none
class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), intent(in) :: form !! Input format code ("XV" or "EL")
integer(I4B), intent(out) :: ierr !! Error code

ierr = 0
end subroutine symba_io_read_frame_info
!module subroutine symba_io_read_frame_info(self, iu, param, form, ierr)
! !! author: David A. Minton
! !!
! !! Reads a single frame of a particle info data from a file.
! implicit none
! class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object
! integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
! class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! character(*), intent(in) :: form !! Input format code ("XV" or "EL")
! integer(I4B), intent(out) :: ierr !! Error code
!
! ierr = 0
! end subroutine symba_io_read_frame_info


module subroutine symba_io_write_discard(self, param)
Expand Down Expand Up @@ -280,13 +346,5 @@ module subroutine symba_io_write_discard(self, param)
return
end subroutine symba_io_write_discard


module subroutine symba_io_write_frame_info(self, iu, param)
implicit none
class(symba_particle_info), intent(in) :: self !! SyMBA particle info object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine symba_io_write_frame_info

end submodule s_symba_io

2 changes: 2 additions & 0 deletions src/symba/symba_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,10 +162,12 @@ module subroutine symba_setup_tp(self, n, param)
if (allocated(self%nplenc)) deallocate(self%nplenc)
if (allocated(self%levelg)) deallocate(self%levelg)
if (allocated(self%levelm)) deallocate(self%levelm)
if (allocated(self%info)) deallocate(self%info)

allocate(self%nplenc(n))
allocate(self%levelg(n))
allocate(self%levelm(n))
allocate(self%info(n))

self%nplenc(:) = 0
self%levelg(:) = -1
Expand Down
Loading

0 comments on commit 88ba918

Please sign in to comment.