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

Commit

Permalink
Browse files Browse the repository at this point in the history 
Improved file i/o robustness and dump parameter file formatting. Now we pass random seeds to the dump param file
  • Loading branch information
@daminton
daminton committed May 25, 2021
1 parent 72d758d commit 0e62487
Show file tree
Hide file tree
Showing 5 changed files with 87 additions and 79 deletions.
13 changes: 6 additions & 7 deletions src/swiftest/swiftest_read_pl_in.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,6 @@
else
self%rhill(:) = 0.0_DP
end if
self%status(:) = ACTIVE
if (param%lclose) then
read(LUN, iostat = ierr) self%radius(:)
else
Expand All @@ -105,21 +104,21 @@
read(LUN, iostat = ierr) self%Ip(:,:)
read(LUN, iostat = ierr) self%rot(:,:)
end if
if (ierr /= 0) then
write(*,*) 'An error occurred reading in ',trim(adjustl(param%inplfile))
call util_exit(FAILURE)
end if
self%status(:) = ACTIVE
end if
close(unit = LUN)
if (ierr /= 0 ) then
write(*,*) 'Error reading in massive body initial conditions from ',trim(adjustl(param%inplfile))
call util_exit(FAILURE)
end if

do i = 1, npl
self%info(1)%origin_type = "Central body"
do i = 2, npl
self%info(i)%origin_xh(:) = self%xh(:,i)
self%info(i)%origin_vh(:) = self%vh(:,i)
self%info(i)%origin_time = 0.0_DP
self%info(i)%origin_type = "Initial conditions"
end do
self%info(1)%origin_type = "Central body"

! Give massive bodies a positive id
self%id(:) = abs(self%id(:))
Expand Down
5 changes: 3 additions & 2 deletions src/swiftest/swiftest_read_tp_in.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,11 @@
write(*,*) 'Error opening test particle initial conditions file ',trim(adjustl(param%intpfile))
return
end if
ntp = 0
if (is_ascii) then
read(lun, *) ntp
read(lun, *, iostat=ierr) ntp
else
read(lun) ntp
read(lun, iostat=ierr) ntp
end if
if (ntp <= 0) return

Expand Down
1 change: 1 addition & 0 deletions src/user/user_dump_param.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ module subroutine user_dump_param(param,t)
param_dump%t0 = t
param_dump%inplfile = trim(adjustl(DUMP_PL_FILE(idx)))
param_dump%intpfile = trim(adjustl(DUMP_TP_FILE(idx)))
param_dump%in_type = "REAL8"
open(unit = LUN, file = DUMP_PARAM_FILE(idx), status='replace', form = 'formatted', iostat =ierr)
if (ierr /=0) then
write(*,*) 'Swiftest error.'
Expand Down
2 changes: 1 addition & 1 deletion src/user/user_read_param_in.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ module subroutine user_read_param_in(param,inparfile)
write(*,*) "EXTRA_FORCE = ",param%lextra_force
write(*,*) "BIG_DISCARD = ",param%lbig_discard
write(*,*) "RHILL_PRESENT = ",param%lrhill_present
write(*,*) "SEED = ",size(param%seed), param%seed(:)
write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:)
ierr = 0

! Added by D. Minton
Expand Down
145 changes: 76 additions & 69 deletions src/user/user_udio_writer.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,88 +19,95 @@ module subroutine user_udio_writer(param, unit, iotype, v_list, iostat, iomsg)
character(len=*), intent(inout) :: iomsg

! Internals
!character(*),parameter :: Ifmt = '(A20,1X,I0/)' !! Format label for integer values
!character(*),parameter :: Rfmt = '(A20,1X,ES25.17/)' !! Format label for real values
!character(*),parameter :: R2fmt = '(A20,2(1X,ES25.17)/)' !! Format label for 2x real values
!character(*),parameter :: Sfmt = '(A20,1X,A/)' !! Format label for string values
!character(*),parameter :: Lfmt = '(A20,1X,L1/)' !! Format label for logical values
!character(*),parameter :: Pfmt = '(A20/)' !! Format label for single parameter string
character(*),parameter :: Ifmt = '(A20,1X,I0)' !! Format label for integer values
character(*),parameter :: Rfmt = '(A20,1X,ES25.17)' !! Format label for real values
character(*),parameter :: R2fmt = '(A20,2(1X,ES25.17))' !! Format label for 2x real values
character(*),parameter :: Sfmt = '(A20,1X,A)' !! Format label for string values
character(*),parameter :: Lfmt = '(A20,1X,L1)' !! Format label for logical values
character(*),parameter :: Pfmt = '(A20)' !! Format label for single parameter string
character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values
character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values
character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values
character(len=*), parameter :: Afmt = '(A25,1X,64(:,A25,1X))'
character(25) :: param_name, param_value
type character_array
character(25) :: value
end type character_array
type(character_array), dimension(:), allocatable :: param_array
integer(I4B) :: i

write(unit, Rfmt) "T0", param%t0
write(unit, Rfmt) "TSTOP",param%tstop
write(unit, Rfmt) "DT",param%dt
write(unit, Sfmt) "PL_IN",trim(adjustl(param%inplfile))
write(unit, Sfmt) "TP_IN",trim(adjustl(param%intpfile))
write(unit, Sfmt) "IN_TYPE",trim(adjustl(param%out_type))
write(param_name, Afmt) "T0"; write(param_value,Rfmt) param%t0; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "TSTOP"; write(param_value, Rfmt) param%tstop; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "DT"; write(param_value, Rfmt) param%dt; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "PL_IN"; write(param_value, Afmt) trim(adjustl(param%inplfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "TP_in"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "IN_TYPE"; write(param_value, Afmt) trim(adjustl(param%in_type)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
if (param%istep_out > 0) then
write(unit, Ifmt) "ISTEP_OUT",param%istep_out
write(unit, Sfmt) "BIN_OUT",trim(adjustl(param%outfile))
write(unit, Sfmt) "PARTICLE_FILE",trim(adjustl(param%particle_file))
write(unit, Sfmt) "OUT_TYPE",trim(adjustl(param%out_type))
write(unit, Sfmt) "OUT_FORM",trim(adjustl(param%out_form))
write(unit, Sfmt) "OUT_STAT","APPEND"
else
write(unit, Pfmt) "!ISTEP_OUT "
write(unit, Pfmt) "!BIN_OUT"
write(unit, Pfmt) "!OUT_TYPE"
write(unit, Pfmt) "!OUT_FORM"
write(unit, Pfmt) "!OUT_STAT"
write(param_name, Afmt) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "BIN_OUT"; write(param_value, Afmt) trim(adjustl(param%outfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "OUT_TYPE"; write(param_value, Afmt) trim(adjustl(param%out_type)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "OUT_FORM"; write(param_value, Afmt) trim(adjustl(param%out_form)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "OUT_STAT"; write(param_value, Afmt) "APPEND"; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
end if
write(unit, Sfmt) "ENC_OUT",trim(adjustl(param%encounter_file))
write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%encounter_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
if (param%istep_dump > 0) then
write(unit, Ifmt) "ISTEP_DUMP",param%istep_dump
else
write(unit, Pfmt) "!ISTEP_DUMP"
write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
end if
if (param%j2rp2 > VSMALL) then
write(unit, Rfmt) "J2 ",param%j2rp2
write(param_name, Afmt) "J2 "; write(param_value, Rfmt) param%j2rp2; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
if (param%j4rp4 > VSMALL) then
write(unit, Rfmt) "J4 ",param%j4rp4
else
write(unit, Pfmt) "!J4 "
write(param_name, Afmt) "J4 "; write(param_value, Rfmt) param%j4rp4; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
end if
else
write(unit, Pfmt) "!J2 "
write(unit, Pfmt) "!J4 "
end if
write(unit, Rfmt) "CHK_RMIN",param%rmin
write(unit, Rfmt) "CHK_RMAX",param%rmax
write(unit, Rfmt) "CHK_EJECT",param%rmaxu
write(unit, Rfmt) "CHK_QMIN",param%qmin
write(param_name, Afmt) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "CHK_EJECT"; write(param_value, Rfmt) param%rmaxu; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "CHK_QMIN"; write(param_value, Rfmt) param%qmin; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
if (param%qmin >= 0.0_DP) then
write(unit, Sfmt) "CHK_QMIN_COORD",trim(adjustl(param%qmin_coord))
write(unit, R2fmt) "CHK_QMIN_RANGE",param%qmin_alo, param%qmin_ahi
else
write(unit, Pfmt) "!CHK_QMIN_COORD"
write(unit, Pfmt) "!CHK_QMIN_RANGE"
write(param_name, Afmt) "CHK_QMIN_COORD"; write(param_value, Afmt) trim(adjustl(param%qmin_coord)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
allocate(param_array(2))
write(param_array(1)%value, Rfmt) param%qmin_alo
write(param_array(2)%value, Rfmt) param%qmin_ahi
write(param_name, Afmt) "CHK_QMIN_RANGE"; write(unit, Afmt) adjustl(param_name), adjustl(param_array(1)%value), adjustl(param_array(2)%value)
end if
if (param%lmtiny) write(unit, Rfmt) "MTINY",param%mtiny
write(unit, Rfmt) "MU2KG",MU2KG
write(unit, Rfmt) "TU2S",TU2S
write(unit, Rfmt) "DU2M",DU2M
write(unit,*) "SEED",size(param%seed),param%seed(:)

write(unit, Lfmt) "EXTRA_FORCE",param%lextra_force
write(unit, Lfmt) "BIG_DISCARD",param%lbig_discard
write(unit, Lfmt) "RHILL_PRESENT",param%lrhill_present
write(unit, Lfmt) "CHK_CLOSE",param%lclose
write(unit, Lfmt) "FRAGMENTATION", param%lfragmentation
write(unit, Lfmt) "ENERGY", param%lenergy
write(unit, Lfmt) "ROTATION", param%lrotation
write(unit, Lfmt) "TIDES", param%ltides
write(unit, Lfmt) "GR", param%lgr
write(unit, Lfmt) "YARKOVSKY", param%lyarkovsky
write(unit, Lfmt) "YORP", param%lyorp
write(unit, Lfmt) "RINGMOONS", param%lringmoons
if (param%lringmoons) write(unit, Sfmt) "RING_OUTFILE",trim(adjustl(param%ring_outfile))
if (param%lmtiny) then
write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
end if
write(param_name, Afmt) "MU2KG"; write(param_value, Rfmt) MU2KG; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "TU2S"; write(param_value, Rfmt) TU2S ; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "DU2M"; write(param_value, Rfmt) DU2M; write(unit, Afmt) adjustl(param_name), adjustl(param_value)

! Special handling is required for writing the random number seed array as its size is not known until runtime
! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements
write(param_name, Afmt) "SEED"
if (allocated(param_array)) deallocate(param_array)
allocate(param_array(0:size(param%seed)))
write(param_array(0)%value, Ifmt) size(param%seed)
do i = 1, size(param%seed)
write(param_array(i)%value, Ifmt) param%seed(i)
end do
write(unit, Afmt, advance='no') adjustl(param_name), adjustl(param_array(0)%value)
do i = 1, size(param%seed)
if (i < size(param%seed)) then
write(unit, Afmt, advance='no') adjustl(param_array(i)%value)
else
write(unit, Afmt) adjustl(param_array(i)%value)
end if
end do
! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

write(param_name, Afmt) "EXTRA_FORCE"; write(param_value, Lfmt) param%lextra_force; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "BIG_DISCARD"; write(param_value, Lfmt) param%lbig_discard; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "RHILL_PRESENT"; write(param_value, Lfmt) param%lrhill_present; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "CHK_CLOSE"; write(param_value, Lfmt) param%lclose; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "ENERGY"; write(param_value, Lfmt) param%lenergy; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "GR"; write(param_value, Lfmt) param%lgr; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "YARKOVSKY"; write(param_value, Lfmt) param%lyarkovsky; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "YORP"; write(param_value, Lfmt) param%lyorp; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "RINGMOONS"; write(param_value, Lfmt) param%lringmoons; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
if (param%lringmoons) then
write(param_name, Afmt) "RING_OUTFILE"; write(param_value, Afmt) trim(adjustl(param%ring_outfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
end if

iostat = 0

return

Expand Down

0 comments on commit 0e62487

Please sign in to comment.