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 OOPTides
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 19, 2021
2 parents a841ad0 + 7182995 commit 7c92ca2
Show file tree
Hide file tree
Showing 10 changed files with 324 additions and 192 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
!README.md
!README.swifter
!*.in
dump*
!example/cleanup
!example/swifter_symba_omp
!Makefile
Expand Down
381 changes: 249 additions & 132 deletions src/io/io.f90
View file

Large diffs are not rendered by default.

16 changes: 12 additions & 4 deletions src/main/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ program swiftest_driver
integer(I8B) :: iout !! Output cadence counter
integer(I8B) :: nloops !! Number of steps to take in the simulation
integer(I4B) :: iu !! Unit number of binary file
real(DP) :: old_t_final = 0.0_DP !! Output time at which writing should start, in order to prevent duplicate lines being written for restarts

ierr = io_get_args(integrator, param_file_name)
if (ierr /= 0) then
Expand All @@ -35,7 +36,7 @@ program swiftest_driver
param%integrator = integrator

call setup_construct_system(nbody_system, param)
call param%read_from_file(param_file_name)
call param%read_in(param_file_name)

associate(t => param%t, &
t0 => param%t0, &
Expand All @@ -50,8 +51,15 @@ program swiftest_driver
iout = istep_out
idump = istep_dump
nloops = ceiling(tstop / dt, kind=I8B)
if (istep_out > 0) call nbody_system%write_frame(iu, param)
call nbody_system%dump(param)
! Prevent duplicate frames from being written if this is a restarted run
if (param%lrestart) then
old_t_final = nbody_system%get_old_t_final(param)
else
old_t_final = t0
if (istep_out > 0) call nbody_system%write_frame(iu, param)
call nbody_system%dump(param)
end if


!> Define the maximum number of threads
nthreads = 1 ! In the *serial* case
Expand All @@ -73,7 +81,7 @@ program swiftest_driver
if (istep_out > 0) then
iout = iout - 1
if (iout == 0) then
call nbody_system%write_frame(iu, param)
if (t > old_t_final) call nbody_system%write_frame(iu, param)
iout = istep_out
end if
end if
Expand Down
78 changes: 40 additions & 38 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,8 @@ module swiftest_classes
character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body
character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies
character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles
character(STRMAX) :: in_type = ASCII_TYPE !! Format of input data files
character(STRMAX) :: in_type = ASCII_TYPE !! Data representation type of input data files
character(STRMAX) :: in_form = XV !! Format of input data files (EL or XV)
integer(I4B) :: istep_out = -1 !! Number of time steps between binary outputs
character(STRMAX) :: outfile = BIN_OUTFILE !! Name of output binary file
character(STRMAX) :: out_type = REAL8_TYPE !! Binary format of output file
Expand Down Expand Up @@ -78,10 +79,10 @@ module swiftest_classes
logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect
logical :: lyorp = .false. !! Turn on YORP effect
contains
procedure :: reader => io_param_reader
procedure :: writer => io_param_writer
procedure :: dump => io_dump_param
procedure :: read_from_file => io_read_param_in
procedure :: reader => io_param_reader
procedure :: writer => io_param_writer
procedure :: dump => io_dump_param
procedure :: read_in => io_read_in_param
end type swiftest_parameters

!********************************************************************************************************************************
Expand All @@ -93,7 +94,6 @@ module swiftest_classes
contains
!! The minimal methods that all systems must have
procedure :: dump => io_dump_swiftest
procedure(abstract_initialize), deferred :: initialize
procedure(abstract_read_frame), deferred :: read_frame
procedure(abstract_write_frame), deferred :: write_frame
end type swiftest_base
Expand Down Expand Up @@ -128,7 +128,7 @@ module swiftest_classes
real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body
real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body
contains
procedure :: initialize => io_read_cb_in !! I/O routine for reading in central body data
procedure :: read_in => io_read_in_cb !! I/O routine for reading in central body data
procedure :: read_frame => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body
procedure :: write_frame => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body
end type swiftest_cb
Expand Down Expand Up @@ -174,7 +174,7 @@ module swiftest_classes
procedure :: drift => drift_body !! Loop through bodies and call Danby drift routine on heliocentric variables
procedure :: v2pv => gr_vh2pv_body !! Converts from velocity to psudeovelocity for GR calculations using symplectic integrators
procedure :: pv2v => gr_pv2vh_body !! Converts from psudeovelocity to velocity for GR calculations using symplectic integrators
procedure :: initialize => io_read_body_in !! Read in body initial conditions from a file
procedure :: read_in => io_read_in_body !! Read in body initial conditions from a file
procedure :: read_frame => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
procedure :: write_frame => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
procedure :: accel_obl => obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
Expand Down Expand Up @@ -306,6 +306,7 @@ module swiftest_classes
procedure :: discard => discard_system !! Perform a discard step on the system
procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time
procedure :: dump => io_dump_system !! Dump the state of the system to a file
procedure :: get_old_t_final => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output.
procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file
procedure :: write_discard => io_write_discard !! Write out information about discarded test particles
procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file
Expand Down Expand Up @@ -356,12 +357,6 @@ subroutine abstract_accel(self, system, param, t, lbeg)
logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step
end subroutine abstract_accel

subroutine abstract_initialize(self, param)
import swiftest_base, swiftest_parameters
class(swiftest_base), intent(inout) :: self !! Swiftest base object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine abstract_initialize

subroutine abstract_kick_body(self, system, param, t, dt, lbeg)
import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP
implicit none
Expand All @@ -373,13 +368,13 @@ subroutine abstract_kick_body(self, system, param, t, dt, lbeg)
logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not.
end subroutine abstract_kick_body

subroutine abstract_read_frame(self, iu, param, form)
function abstract_read_frame(self, iu, param) result(ierr)
import DP, I4B, swiftest_base, swiftest_parameters
class(swiftest_base), intent(inout) :: self !! Swiftest base object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), intent(in) :: form !! Input format code ("XV" or "EL")
end subroutine abstract_read_frame
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful
end function abstract_read_frame

subroutine abstract_set_mu(self, cb)
import swiftest_body, swiftest_cb
Expand Down Expand Up @@ -579,6 +574,13 @@ module function io_get_args(integrator, param_file_name) result(ierr)
integer(I4B) :: ierr !! I/O error code
end function io_get_args

module function io_get_old_t_final_system(self, param) result(old_t_final)
implicit none
class(swiftest_nbody_system), intent(in) :: self
class(swiftest_parameters), intent(in) :: param
real(DP) :: old_t_final
end function io_get_old_t_final_system

module function io_get_token(buffer, ifirst, ilast, ierr) result(token)
implicit none
character(len=*), intent(in) :: buffer !! Input string buffer
Expand Down Expand Up @@ -610,47 +612,47 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0
end subroutine io_param_writer

module subroutine io_read_body_in(self, param)
module subroutine io_read_in_body(self, param)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine io_read_body_in
end subroutine io_read_in_body

module subroutine io_read_cb_in(self, param)
module subroutine io_read_in_cb(self, param)
implicit none
class(swiftest_cb), intent(inout) :: self !! Swiftest central body object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine io_read_cb_in
end subroutine io_read_in_cb

module subroutine io_read_param_in(self, param_file_name)
module subroutine io_read_in_param(self, param_file_name)
implicit none
class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
end subroutine io_read_param_in
end subroutine io_read_in_param

module subroutine io_read_frame_body(self, iu, param, form)
module function io_read_frame_body(self, iu, param) result(ierr)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), intent(in) :: form !! Input format code ("XV" or "EL")
end subroutine io_read_frame_body
class(swiftest_body), intent(inout) :: self !! Swiftest body object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful
end function io_read_frame_body

module subroutine io_read_frame_cb(self, iu, param, form)
module function io_read_frame_cb(self, iu, param) result(ierr)
implicit none
class(swiftest_cb), intent(inout) :: self !! Swiftest central body object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), intent(in) :: form !! Input format code ("XV" or "EL")
end subroutine io_read_frame_cb
class(swiftest_cb), intent(inout) :: self !! Swiftest central body object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful
end function io_read_frame_cb

module subroutine io_read_frame_system(self, iu, param, form)
module function io_read_frame_system(self, iu, param) result(ierr)
implicit none
class(swiftest_nbody_system),intent(inout) :: self !! Swiftest system object
integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
character(*), intent(in) :: form !! Input format code ("XV" or "EL")
end subroutine io_read_frame_system
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful
end function io_read_frame_system

module subroutine io_write_discard(self, param)
implicit none
Expand Down
2 changes: 1 addition & 1 deletion src/modules/swiftest_globals.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ module swiftest_globals
character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ]
character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ]
character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ]
character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.dat', 'dump_param2.dat']
character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.in', 'dump_param2.in']

!> Default file names that can be changed by the user in the parameters file
character(*), parameter :: CB_INFILE = 'cb.in'
Expand Down
10 changes: 6 additions & 4 deletions src/setup/setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,16 +133,18 @@ module subroutine setup_initialize_system(self, param)
! Arguments
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
call self%cb%initialize(param)
call self%pl%initialize(param)
call self%tp%initialize(param)

call self%cb%read_in(param)
call self%pl%read_in(param)
call self%tp%read_in(param)
call self%validate_ids(param)
call self%set_msys()
call self%pl%set_mu(self%cb)
call self%tp%set_mu(self%cb)
call self%pl%eucl_index()
if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb)
self%pl%lfirst = param%lfirstkick
self%tp%lfirst = param%lfirstkick
return
end subroutine setup_initialize_system

Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body)
system%Mescape = system%Mescape + pl%mass(ipl)
do i = 1, pl%nbody
if (i == ipl) cycle
pe = pe - pl%mass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i))
pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i))
end do

Ltot(:) = 0.0_DP
Expand Down
3 changes: 2 additions & 1 deletion src/symba/symba_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ module subroutine symba_io_read_particle(system, param)
associate(npl => pl%nbody, ntp => tp%nbody)
do
lmatch = .false.
read(LUN, err = 667, iomsg = errmsg) id
read(LUN, err = 667, iomsg = errmsg, end = 333) id

if (idx == cb%id) then
read(LUN, err = 667, iomsg = errmsg) cb%info
Expand Down Expand Up @@ -308,6 +308,7 @@ module subroutine symba_io_read_particle(system, param)
end select
end select

333 continue
return

667 continue
Expand Down
7 changes: 4 additions & 3 deletions src/symba/symba_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,8 @@ module subroutine symba_step_system(self, param, t, dt)
call self%reset(param)
lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt, 0)
if (lencounter) then
tp%lfirst = pl%lfirst
call self%interp(param, t, dt)
pl%lfirst = .true.
tp%lfirst = .true.
param%lfirstkick = .true.
else
self%irec = -1
call helio_step_system(self, param, t, dt)
Expand Down Expand Up @@ -266,6 +264,9 @@ module subroutine symba_step_reset_system(self, param)

call system%pl_adds%setup(0, param)
call system%pl_discards%setup(0, param)

tp%lfirst = param%lfirstkick
pl%lfirst = param%lfirstkick
end associate
end select
end select
Expand Down
16 changes: 8 additions & 8 deletions src/symba/symba_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -455,15 +455,15 @@ module subroutine symba_util_rearray_pl(self, system, param)
end if

! Reset all of the status flags for this body
where(pl%status(:) /= INACTIVE)
pl%status(:) = ACTIVE
pl%ldiscard(:) = .false.
pl%lcollision(:) = .false.
pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
pl%lmask(:) = .true.
where(pl%status(1:pl%nbody) /= INACTIVE)
pl%status(1:pl%nbody) = ACTIVE
pl%ldiscard(1:pl%nbody) = .false.
pl%lcollision(1:pl%nbody) = .false.
pl%lmtiny(1:pl%nbody) = pl%Gmass(1:pl%nbody) > param%GMTINY
pl%lmask(1:pl%nbody) = .true.
elsewhere
pl%ldiscard(:) = .true.
pl%lmask(:) = .false.
pl%ldiscard(1:pl%nbody) = .true.
pl%lmask(1:pl%nbody) = .false.
end where

pl%nplm = count(pl%lmtiny(1:pl%nbody) .and. pl%lmask(1:pl%nbody))
Expand Down

0 comments on commit 7c92ca2

Please sign in to comment.