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

Commit

Permalink
Merge branch 'debug' into FastEncounterCheck
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Sep 27, 2021
2 parents bff585f + 968472a commit f97ca1c
Show file tree
Hide file tree
Showing 10 changed files with 534 additions and 120 deletions.
2 changes: 1 addition & 1 deletion examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@
print(f'TP_IN {swifter_tp}')
print(f'IN_TYPE ASCII')
print(f'ISTEP_OUT {iout:d}')
print(f'ISTEP_DUMP {iout:d}')
print(f'ISTEP_DUMP {10*iout:d}')
print(f'BIN_OUT {swifter_bin}')
print(f'OUT_TYPE REAL8')
print(f'OUT_FORM XV')
Expand Down
2 changes: 1 addition & 1 deletion examples/rmvs_swifter_comparison/1pl_1tp_encounter/param.swiftest.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ PL_IN pl.swiftest.in
TP_IN tp.swiftest.in
IN_TYPE REAL8
ISTEP_OUT 1
ISTEP_DUMP 1
ISTEP_DUMP 10
BIN_OUT bin.swiftest.nc
OUT_TYPE NETCDF_DOUBLE
OUT_FORM XVEL
Expand Down
16 changes: 8 additions & 8 deletions examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb
View file

Large diffs are not rendered by default.

163 changes: 110 additions & 53 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,18 @@ module subroutine io_conservation_report(self, param, lterminal)
"; DM/M0 = ", ES12.5)'

associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody)
if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE) .and. (param%energy_out /= "")) then
if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then
if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then
open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write", err = 667, iomsg = errmsg)
write(EGYIU,EGYHEADER, err = 667, iomsg = errmsg)
else
open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg)
end if
end if

call pl%vb2vh(cb)
call pl%xh2xb(cb)

call system%get_energy_and_momentum(param)
ke_orbit_now = system%ke_orbit
ke_spin_now = system%ke_spin
Expand All @@ -56,7 +58,7 @@ module subroutine io_conservation_report(self, param, lterminal)
param%lfirstenergy = .false.
end if

if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE) .and. (param%energy_out /= "")) then
if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then
write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, param%Ecollisions, Ltot_now, GMtot_now
close(EGYIU, err = 667, iomsg = errmsg)
end if
Expand All @@ -71,11 +73,11 @@ module subroutine io_conservation_report(self, param, lterminal)
if (abs(Merror) > 100 * epsilon(Merror)) then
write(*,*) "Severe error! Mass not conserved! Halting!"
call pl%xv2el(cb)
call param%nciu%open(param)
!call param%nciu%open(param)
call self%write_hdr(param%nciu, param)
call cb%write_frame(param%nciu, param)
call pl%write_frame(param%nciu, param)
call param%nciu%close(param)
call param%nciu%close()
call util_exit(FAILURE)
end if
end if
Expand Down Expand Up @@ -200,9 +202,7 @@ module subroutine io_dump_particle_info_base(self, param, idx)
close(unit = LUN, err = 667, iomsg = errmsg)
!else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
call param%nciu%open(param)
call self%write_particle_info(param%nciu)
call param%nciu%close(param)
end if

return
Expand Down Expand Up @@ -275,41 +275,52 @@ module subroutine io_dump_system(self, param)
real(DP) :: tfrac
character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active pl, tp = ", I5, ", ", I5)'
character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)'
logical, save :: lfirst = .true.

if (lfirst) then
lfirst = .false.
if (param%lenergy) call self%conservation_report(param, lterminal=.false.)
else
allocate(dump_param, source=param)
param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx)))

allocate(dump_param, source=param)
param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx)))
dump_param%in_form = XV
dump_param%out_form = XV
dump_param%out_stat = 'APPEND'
if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then
dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx)))
dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx)))
dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx)))
dump_param%in_form = XV
dump_param%out_stat = 'APPEND'
dump_param%in_type = REAL8_TYPE
dump_param%T0 = param%t
else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx)))
dump_param%in_type = NETCDF_DOUBLE_TYPE
end if
dump_param%T0 = param%t
dump_param%ioutput = 0

call dump_param%dump(param_file_name)
call dump_param%dump(param_file_name)

dump_param%out_form = XV
if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then
call self%cb%dump(dump_param)
call self%pl%dump(dump_param)
call self%tp%dump(dump_param)
else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
call dump_param%nciu%initialize(dump_param)
call self%write_hdr(dump_param%nciu, dump_param)
call self%cb%write_frame(dump_param%nciu, dump_param)
call self%pl%write_frame(dump_param%nciu, dump_param)
call self%tp%write_frame(dump_param%nciu, dump_param)
call dump_param%nciu%close()
! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk)
call param%nciu%sync()
end if

idx = idx + 1
if (idx > NDUMPFILES) idx = 1
idx = idx + 1
if (idx > NDUMPFILES) idx = 1

tfrac = (param%t - param%t0) / (param%tstop - param%t0)

select type(pl => self%pl)
class is (symba_pl)
write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, self%tp%nbody
class default
write(*, statusfmt) param%t, tfrac, pl%nbody, self%tp%nbody
end select
end if
tfrac = (param%t - param%t0) / (param%tstop - param%t0)

select type(pl => self%pl)
class is (symba_pl)
write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, self%tp%nbody
class default
write(*, statusfmt) param%t, tfrac, pl%nbody, self%tp%nbody
end select

return
end subroutine io_dump_system
Expand Down Expand Up @@ -393,7 +404,7 @@ module function io_get_old_t_final_system(self, param) result(old_t_final)
class(swiftest_nbody_system), allocatable :: tmpsys
class(swiftest_parameters), allocatable :: tmpparam
integer(I4B) :: ierr, iu = LUN
character(len=STRMAX) :: errmsg
character(len=STRMAX) :: errmsg

old_t_final = 0.0_DP
allocate(tmpsys, source=self)
Expand Down Expand Up @@ -1143,7 +1154,42 @@ module subroutine io_param_writer_one_QP(param_name, param_value, unit)
end subroutine io_param_writer_one_QP


module subroutine io_read_in_body(self, param)
module subroutine io_read_in_base(self,param)
!! author: Carlisle A. Wishard and David A. Minton
!!
!! Reads in either a central body, test particle, or massive body object. For the swiftest_body types (non-central body), it allocates array space for them
implicit none
class(swiftest_base), intent(inout) :: self !! Swiftest base object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful

select case(param%in_type)
case(NETCDF_DOUBLE_TYPE, NETCDF_FLOAT_TYPE)
select type(self)
class is (swiftest_body)
if (self%nbody == 0) return
call self%setup(self%nbody, param)
end select

ierr = self%read_frame(param%nciu, param)
if (ierr == 0) return
case default
select type(self)
class is (swiftest_body)
call io_read_in_body(self, param)
class is (swiftest_cb)
call io_read_in_cb(self, param)
end select
return
end select

667 continue
write(*,*) "Error reading body in io_read_in_base"
end subroutine io_read_in_base


subroutine io_read_in_body(self, param)
!! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
!!
!! Read in either test particle or massive body data
Expand All @@ -1157,20 +1203,19 @@ module subroutine io_read_in_body(self, param)
! Internals
integer(I4B) :: iu = LUN
integer(I4B) :: i, nbody
logical :: is_ascii, is_pl
logical :: is_ascii
character(len=:), allocatable :: infile
real(DP) :: t
character(STRMAX) :: errmsg
integer(I4B) :: ierr

! Select the appropriate polymorphic class (test particle or massive body)

select type(self)
class is (swiftest_pl)
infile = param%inplfile
is_pl = .true.
class is (swiftest_tp)
infile = param%intpfile
is_pl = .false.
end select

is_ascii = (param%in_type == 'ASCII')
Expand Down Expand Up @@ -1207,7 +1252,7 @@ module subroutine io_read_in_body(self, param)
end subroutine io_read_in_body


module subroutine io_read_in_cb(self, param)
subroutine io_read_in_cb(self, param)
!! author: David A. Minton
!!
!! Reads in central body data
Expand Down Expand Up @@ -1268,6 +1313,30 @@ module subroutine io_read_in_cb(self, param)
end subroutine io_read_in_cb


module subroutine io_read_in_system(self, param)
!! author: David A. Minton and Carlisle A. Wishard
!!
!! Reads in the system from input files
implicit none
! Arguments
class(swiftest_nbody_system), intent(inout) :: self
class(swiftest_parameters), intent(inout) :: param
! Internals
integer(I4B) :: ierr

if ((param%in_type == NETCDF_DOUBLE_TYPE) .or. (param%in_type == NETCDF_FLOAT_TYPE)) then
ierr = self%read_frame(param%nciu, param)
if (ierr /=0) call util_exit(FAILURE)
else
call self%cb%read_in(param)
call self%pl%read_in(param)
call self%tp%read_in(param)
end if

return
end subroutine io_read_in_system


function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, &
xh1, xh2, vh1, vh2, enc_out, out_type) result(ierr)
!! author: David A. Minton
Expand Down Expand Up @@ -1787,9 +1856,7 @@ module subroutine io_write_discard(self, param)

! Record the discarded body metadata information to file
if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
call param%nciu%open(param)
call tp_discards%write_particle_info(param%nciu)
call param%nciu%close(param)
end if

if (param%discard_out == "") return
Expand Down Expand Up @@ -2094,29 +2161,20 @@ module subroutine io_write_frame_system(self, param)
errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN"
goto 667
end if
call param%nciu%open(param)
case('NEW')
if (fileExists) then
errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = OLD, REPLACE, or UNKNOWN"
errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN"
goto 667
end if
call param%nciu%initialize(param)
case('REPLACE', 'UNKNOWN')
if (fileExists) then
open(file=param%outfile, unit=iu, status='OLD')
close (unit=BINUNIT, status="delete")
end if
end select

select case(param%out_stat)
case('APPEND')
call param%nciu%open(param)
case('NEW', 'REPLACE', 'UNKNOWN')
call param%nciu%initialize(param)
call param%nciu%close(param)
call param%nciu%open(param)
end select

lfirst = .false.
else
call param%nciu%open(param)
!call param%nciu%open(param)
end if
end if

Expand All @@ -2142,7 +2200,6 @@ module subroutine io_write_frame_system(self, param)
call cb%write_frame(param%nciu, param)
call pl%write_frame(param%nciu, param)
call tp%write_frame(param%nciu, param)
call param%nciu%close(param)
end if

return
Expand Down
8 changes: 6 additions & 2 deletions src/main/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,11 +57,13 @@ program swiftest_driver
ioutput = ioutput_t0
! Prevent duplicate frames from being written if this is a restarted run
if ((param%lrestart) .and. ((param%out_type == REAL8_TYPE) .or. param%out_type == REAL4_TYPE)) then
old_t_final = nbody_system%get_old_t_final(param)
old_t_final = nbody_system%get_old_t_final_bin(param)
else if ((param%lrestart) .and. ((param%out_type == NETCDF_DOUBLE_TYPE) .or. param%out_type == NETCDF_FLOAT_TYPE)) then
old_t_final = nbody_system%get_old_t_final_netcdf(param)
else
old_t_final = t0
if (istep_out > 0) call nbody_system%write_frame(param)
call nbody_system%dump(param)
if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum
end if

!> Define the maximum number of threads
Expand Down Expand Up @@ -104,6 +106,8 @@ program swiftest_driver
end do
end associate

call nbody_system%finalize(param)

call util_exit(SUCCESS)

stop
Expand Down
Loading

0 comments on commit f97ca1c

Please sign in to comment.