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

Commit

Permalink
Merge branch 'coarray_test_particles' of github.itap.purdue.edu:Minto…
Browse files Browse the repository at this point in the history 
…nGroup/swiftest into coarray_test_particles
  • Loading branch information
@daminton
daminton committed Mar 20, 2023
2 parents 5df51a6 + 555c116 commit 79af8c5
Show file tree
Hide file tree
Showing 8 changed files with 92 additions and 106 deletions.
5 changes: 0 additions & 5 deletions cmake/Modules/SetFortranFlags.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,11 +325,6 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_RELEASE}"
"-pg -fbacktrace"
)

SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_PROFILE}"
Fortran "-check bounds,pointers,uninit" # Intel
"-fcheck=bounds,pointer,mem"
)

# Sanitize
SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
Fortran "-fsanitize=address,undefined" # Gnu
Expand Down
8 changes: 4 additions & 4 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -499,7 +499,7 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vp
end subroutine encounter_check_all_triangular_plplm


subroutine encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, &
subroutine encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, &
nenc, index1, index2, lvdotr)
!! author: David A. Minton
!!
Expand All @@ -511,7 +511,7 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, ren
integer(I4B), intent(in) :: ntp !! Total number of test particles
real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies
real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies
real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of massive bodies
real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of massive bodies
real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies
real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter
real(DP), intent(in) :: dt !! Step size
Expand All @@ -529,12 +529,12 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, ren
renct(:) = 0.0_DP

!$omp parallel do default(private) schedule(dynamic)&
!$omp shared(rpl, vpl, xtp, vtp, renc, renct, lenc, ind_arr) &
!$omp shared(rpl, vpl, rtp, vtp, renc, renct, lenc, ind_arr) &
!$omp firstprivate(npl, ntp, dt)
do i = 1, npl
call encounter_check_all_triangular_one(0, ntp, rpl(1,i), rpl(2,i), rpl(3,i), &
vpl(1,i), vpl(2,i), vpl(3,i), &
xtp(1,:), xtp(2,:), xtp(3,:), &
rtp(1,:), rtp(2,:), rtp(3,:), &
vtp(1,:), vtp(2,:), vtp(3,:), &
renc(i), renct(:), dt, ind_arr(:), lenc(i))
if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i
Expand Down
133 changes: 65 additions & 68 deletions src/swiftest/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ program swiftest_driver
use swiftest
implicit none

class(base_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated
class(swiftest_nbody_system), allocatable :: nbody_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 globals for symbolic names)
character(len=:), allocatable :: param_file_name !! Name of the file containing user-defined parameters
Expand Down Expand Up @@ -70,78 +70,75 @@ program swiftest_driver

! Construct the main n-body nbody_system using the user-input integrator to choose the type of nbody_system
call swiftest_util_setup_construct_system(nbody_system, param)
select type(nbody_system)
class is (swiftest_nbody_system)

!> Define the maximum number of threads
nthreads = 1 ! In the *serial* case
!$ nthreads = omp_get_max_threads() ! In the *parallel* case
!$ write(param%display_unit,'(a)') ' OpenMP parameters:'
!$ write(param%display_unit,'(a)') ' ------------------'
!$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads
!$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads

call nbody_system%initialize(param)

associate (system_history => nbody_system%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.
call nbody_system%display_run_information(param, integration_timer, phase="first")
if (param%lenergy) then
if (param%lrestart) then
call nbody_system%get_t0_values(param)
else
call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum
end if
call nbody_system%conservation_report(param, lterminal=.true.)

!> Define the maximum number of threads
nthreads = 1 ! In the *serial* case
!$ nthreads = omp_get_max_threads() ! In the *parallel* case
!$ write(param%display_unit,'(a)') ' OpenMP parameters:'
!$ write(param%display_unit,'(a)') ' ------------------'
!$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads
!$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads

call nbody_system%initialize(param)

associate (system_history => nbody_system%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.
call nbody_system%display_run_information(param, integration_timer, phase="first")
if (param%lenergy) then
if (param%lrestart) then
call nbody_system%get_t0_values(param)
else
call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum
end if
call system_history%take_snapshot(param,nbody_system)
call nbody_system%dump(param)

do iloop = istart, nloops
!> Step the nbody_system forward in time
call integration_timer%start()
call nbody_system%step(param, nbody_system%t, dt)
call integration_timer%stop()

nbody_system%t = t0 + iloop * dt

!> Evaluate any discards or collisional outcomes
call nbody_system%discard(param)

!> 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) .or. (iloop == nloops)) then
iout = 0
idump = idump + 1
if (ltstretch) then
nout = nout + 1
istep = floor(istep_out * fstep_out**nout, kind=I4B)
end if

call system_history%take_snapshot(param,nbody_system)

if (idump == dump_cadence) then
idump = 0
call nbody_system%dump(param)
end if

call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out)
call nbody_system%display_run_information(param, integration_timer)
call integration_timer%reset()
if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)
call nbody_system%conservation_report(param, lterminal=.true.)
end if
call system_history%take_snapshot(param,nbody_system)
call nbody_system%dump(param)

do iloop = istart, nloops
!> Step the nbody_system forward in time
call integration_timer%start()
call nbody_system%step(param, nbody_system%t, dt)
call integration_timer%stop()

nbody_system%t = t0 + iloop * dt

!> Evaluate any discards or collisional outcomes
call nbody_system%discard(param)

!> 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) .or. (iloop == nloops)) then
iout = 0
idump = idump + 1
if (ltstretch) then
nout = nout + 1
istep = floor(istep_out * fstep_out**nout, kind=I4B)
end if

call system_history%take_snapshot(param,nbody_system)

if (idump == dump_cadence) then
idump = 0
call nbody_system%dump(param)
end if

call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out)
call nbody_system%display_run_information(param, integration_timer)
call integration_timer%reset()
if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)

end if
end if

end do
! Dump any remaining history if it exists
call nbody_system%dump(param)
call system_history%dump(param)
call nbody_system%display_run_information(param, integration_timer, phase="last")

end associate
end select
end do
! Dump any remaining history if it exists
call nbody_system%dump(param)
call system_history%dump(param)
call nbody_system%display_run_information(param, integration_timer, phase="last")

end associate
end associate

call base_util_exit(SUCCESS)
Expand Down
2 changes: 0 additions & 2 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1267,8 +1267,6 @@ module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param)
call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, Gmtemp, start=[1,tslot], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar Gmass_varid" )
where(Gmtemp(:) /= Gmtemp(:)) Gmtemp(:) = 0.0_DP
plmmask(:) = plmask(:) .and. Gmtemp(:) > param%GMTINY
else
plmmask(:) = plmask(:)
end if

status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid)
Expand Down
21 changes: 10 additions & 11 deletions src/swiftest/swiftest_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,8 +342,7 @@ module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass
end subroutine swiftest_kick_getacch_int_all_tri_norad_pl



module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, xtp, rpl, GMpl, lmask, acc)
module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, rtp, rpl, GMpl, lmask, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies with parallelisim
Expand All @@ -353,27 +352,27 @@ module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, xtp, rpl, GMpl, lma
implicit none
integer(I4B), intent(in) :: ntp !! Number of test particles
integer(I4B), intent(in) :: npl !! Number of massive bodies
real(DP), dimension(:,:), intent(in) :: xtp !! Test particle position vector array
real(DP), dimension(:,:), intent(in) :: rtp !! Test particle position vector array
real(DP), dimension(:,:), intent(in) :: rpl !! Massive body particle position vector array
real(DP), dimension(:), intent(in) :: GMpl !! Array of massive body G*mass
logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which test particles should be computed
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
! Internals
real(DP) :: rji2
real(DP) :: xr, yr, zr
real(DP) :: rx, ry, rz
integer(I4B) :: i, j

!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, ntp, lmask, xtp, rpl, GMpl) &
!$omp shared(npl, ntp, lmask, rtp, rpl, GMpl) &
!$omp reduction(-:acc)
do i = 1, ntp
if (lmask(i)) then
do j = 1, npl
xr = xtp(1, i) - rpl(1, j)
yr = xtp(2, i) - rpl(2, j)
zr = xtp(3, i) - rpl(3, j)
rji2 = xr**2 + yr**2 + zr**2
call swiftest_kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(j), acc(1,i), acc(2,i), acc(3,i))
do concurrent (j = 1:npl)
rx = rtp(1, i) - rpl(1, j)
ry = rtp(2, i) - rpl(2, j)
rz = rtp(3, i) - rpl(3, j)
rji2 = rx**2 + ry**2 + rz**2
call swiftest_kick_getacch_int_one_tp(rji2, rx, ry, rz, GMpl(j), acc(1,i), acc(2,i), acc(3,i))
end do
end if
end do
Expand Down
4 changes: 2 additions & 2 deletions src/swiftest/swiftest_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -924,11 +924,11 @@ module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine swiftest_kick_getacch_int_all_tri_norad_pl

module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, xtp, rpl, GMpl, lmask, acc)
module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, rtp, rpl, GMpl, lmask, acc)
implicit none
integer(I4B), intent(in) :: ntp !! Number of test particles
integer(I4B), intent(in) :: npl !! Number of massive bodies
real(DP), dimension(:,:), intent(in) :: xtp !! Test particle position vector array
real(DP), dimension(:,:), intent(in) :: rtp !! Test particle position vector array
real(DP), dimension(:,:), intent(in) :: rpl !! Massive body particle position vector array
real(DP), dimension(:), intent(in) :: GMpl !! Array of massive body G*mass
logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which test particles should be computed
Expand Down
21 changes: 9 additions & 12 deletions src/walltime/walltime_implementations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ module subroutine walltime_stop(self)
end subroutine walltime_stop


module subroutine walltime_report(self, message, unit, nsubsteps)
module subroutine walltime_report(self, message, unit)
!! author: David A. Minton
!!
!! Prints the elapsed time information to the terminal
Expand All @@ -47,9 +47,7 @@ module subroutine walltime_report(self, message, unit, nsubsteps)
class(walltimer), intent(inout) :: self !! Walltimer object
character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output
integer(I4B), intent(in) :: unit !! Output file unit for report text to be directed
integer(I4B), optional, intent(in) :: nsubsteps !! Number of substeps used to compute the time per step
! Internals
character(len=*), parameter :: nosubstepfmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5 '
character(len=*), parameter :: substepfmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5, ";' //&
' Interval wall time/step: ", es12.5'
character(len=STRMAX) :: fmt
Expand All @@ -62,17 +60,14 @@ module subroutine walltime_report(self, message, unit, nsubsteps)

call system_clock(count_now)
count_delta_main = count_now - self%count_start_main
count_delta_step = count_now - self%count_start_step
self%wall_main = count_delta_main / (self%count_rate * 1.0_DP)

count_delta_step = self%count_stop_step - self%count_start_step
self%wall_step = count_delta_step / (self%count_rate * 1.0_DP)
if (present(nsubsteps)) then
self%wall_per_substep = self%wall_step / nsubsteps
fmt = '("' // adjustl(message) // '",' // substepfmt // ')'
write(unit,trim(adjustl(fmt))) self%wall_main, self%wall_step, self%wall_per_substep
else
fmt = '("' // adjustl(message) // '",' // nosubstepfmt // ')'
write(unit,trim(adjustl(fmt))) self%wall_main, self%wall_step
end if
self%wall_per_substep = self%wall_step / self%nsubsteps

fmt = '("' // adjustl(message) // '",' // substepfmt // ')'
write(unit,trim(adjustl(fmt))) self%wall_main, self%wall_step, self%wall_per_substep

return
end subroutine walltime_report
Expand All @@ -90,6 +85,7 @@ module subroutine walltime_reset(self)
self%is_paused = .false.
self%wall_step = 0.0_DP
self%wall_per_substep = 0.0_DP
self%nsubsteps = 0

return
end subroutine walltime_reset
Expand Down Expand Up @@ -136,6 +132,7 @@ module subroutine walltime_start(self)
else ! Start a new step timer
call system_clock(self%count_start_step)
end if
self%nsubsteps = self%nsubsteps + 1

return
end subroutine walltime_start
Expand Down
4 changes: 2 additions & 2 deletions src/walltime/walltime_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ module walltime
integer(I8B) :: count_start_step !! Value of the clock ticker at the start of a timed step
integer(I8B) :: count_stop_step !! Value of the clock ticker at the end of a timed step
integer(I8B) :: count_pause !! Value of the clock ticker at the end of a timed step
integer(I4B) :: nsubsteps !! Number of substeps in an interval (number of times the timer is turned off and back on again)
real(DP) :: wall_step !! Value of the step elapsed time
real(DP) :: wall_main !! Value of the main clock elapsed time
real(DP) :: wall_per_substep !! Value of time per substep
Expand All @@ -44,12 +45,11 @@ module walltime


interface
module subroutine walltime_report(self, message, unit, nsubsteps)
module subroutine walltime_report(self, message, unit)
implicit none
class(walltimer), intent(inout) :: self !! Walltimer object
character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output
integer(I4B), intent(in) :: unit !! Output file unit for report text to be directed
integer(I4B), optional, intent(in) :: nsubsteps !! Number of substeps used to compute the time per step
end subroutine walltime_report

module subroutine walltime_reset(self)
Expand Down

0 comments on commit 79af8c5

Please sign in to comment.