From b078c93ff064f91b02c8692b69688617ce654d93 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 10 Mar 2023 10:16:35 -0500 Subject: [PATCH 1/7] Removed bounds checking from profiling mode. Refactored test particle position variable from xtp to rtp --- cmake/Modules/SetFortranFlags.cmake | 5 ----- src/encounter/encounter_check.f90 | 8 ++++---- src/swiftest/swiftest_kick.f90 | 21 ++++++++++----------- src/swiftest/swiftest_module.f90 | 4 ++-- 4 files changed, 16 insertions(+), 22 deletions(-) diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index 891a878f4..97009a25e 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -328,11 +328,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 diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 105aaf7fe..c2013ce6b 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -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 !! @@ -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 @@ -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 diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index 16bb2ac31..54da2f82a 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -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 @@ -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 diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 7913fd1ae..4ac6cb3c3 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -933,11 +933,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 From a28a064db801f86eba8081c30f1c9370daf48dfd Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 10 Mar 2023 15:30:02 -0500 Subject: [PATCH 2/7] Got rid of unnecessary array assignment for when there are not semi-interacting bodies --- src/swiftest/swiftest_io.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 9406e44c2..32f722a7b 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -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) From a312d725470282aa1c96b3e4c2beb8fc0035e982 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 10 Mar 2023 17:57:40 -0500 Subject: [PATCH 3/7] Simplified the wall timer and fixed a problem that was causing it to add in the elapsed time between when the timer was stopped and when the report was caleld to the interval time. --- src/swiftest/swiftest_driver.f90 | 2 +- src/walltime/walltime_implementations.f90 | 20 +++++++++----------- src/walltime/walltime_module.f90 | 4 ++-- 3 files changed, 12 insertions(+), 14 deletions(-) diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 2da264bea..bfa950222 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -124,7 +124,7 @@ program swiftest_driver call nbody_system%dump(param) end if - call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out) + call integration_timer%report(message="Integration steps:", unit=display_unit) call nbody_system%display_run_information(param, integration_timer) call integration_timer%reset() if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) diff --git a/src/walltime/walltime_implementations.f90 b/src/walltime/walltime_implementations.f90 index dcf1a2a60..aaf18f85b 100644 --- a/src/walltime/walltime_implementations.f90 +++ b/src/walltime/walltime_implementations.f90 @@ -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 @@ -47,7 +47,6 @@ 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, ";' //& @@ -62,17 +61,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 @@ -90,6 +86,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 @@ -106,6 +103,7 @@ module subroutine walltime_start_main(self) call system_clock(self%count_start_main, self%count_rate, self%count_max) self%main_is_started = .true. self%wall_main = 0.0_DP + self%nsubsteps = self%nsubsteps + 1 return end subroutine walltime_start_main diff --git a/src/walltime/walltime_module.f90 b/src/walltime/walltime_module.f90 index 6ccce25b1..fa41bae87 100644 --- a/src/walltime/walltime_module.f90 +++ b/src/walltime/walltime_module.f90 @@ -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 @@ -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) From 7cb3c0c843281b342b538f066b776e8a4e5f51c0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 10 Mar 2023 18:00:06 -0500 Subject: [PATCH 4/7] Put the substep counter where it belongs --- src/walltime/walltime_implementations.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/walltime/walltime_implementations.f90 b/src/walltime/walltime_implementations.f90 index aaf18f85b..0a9824d78 100644 --- a/src/walltime/walltime_implementations.f90 +++ b/src/walltime/walltime_implementations.f90 @@ -48,7 +48,6 @@ module subroutine walltime_report(self, message, unit) 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 ! 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 @@ -103,7 +102,6 @@ module subroutine walltime_start_main(self) call system_clock(self%count_start_main, self%count_rate, self%count_max) self%main_is_started = .true. self%wall_main = 0.0_DP - self%nsubsteps = self%nsubsteps + 1 return end subroutine walltime_start_main @@ -134,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 From cbece6cb76628758a4f6cfb109eacec038e075d2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Mar 2023 11:05:07 -0400 Subject: [PATCH 5/7] Updated the sundiver test with more massive body to test for angular momentum conservation --- examples/solar_impact/sundiver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) mode change 100644 => 100755 examples/solar_impact/sundiver.py diff --git a/examples/solar_impact/sundiver.py b/examples/solar_impact/sundiver.py old mode 100644 new mode 100755 index 02e320f69..c279b4ec4 --- a/examples/solar_impact/sundiver.py +++ b/examples/solar_impact/sundiver.py @@ -46,9 +46,9 @@ q = 0.9 * swiftest.RSun * sim.M2DU a = 0.1 e = 1.0 - q / a -M = 1e-4 * swiftest.MEarth * sim.KG2MU +M = 2e0 * swiftest.MEarth * sim.KG2MU R = (3 * M / (4 * np.pi * density)) ** (1.0 / 3.0) -rot = 2 * sim.init_cond.sel(name="Earth")['rot'] +rot = 4 * sim.init_cond.sel(name="Earth")['rot'] sim.add_body(name="Sundiver", a=a, e=e, inc=0.0, capom=0.0, omega=0.0, capm=180.0, mass=M, radius=R, Ip=[0.4,0.4,0.4], rot=rot) sim.get_parameter() From 27392867712d36bdbd7bd256e851c22b65138518 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Mar 2023 11:05:22 -0400 Subject: [PATCH 6/7] Added status check to vb2vh conversion --- src/swiftest/swiftest_util.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index cc1b9e398..fcf130f83 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -494,7 +494,7 @@ module subroutine swiftest_util_coord_vb2vh_pl(self, cb) associate(pl => self, npl => self%nbody) cb%vb(:) = 0.0_DP do i = npl, 1, -1 - cb%vb(:) = cb%vb(:) - pl%Gmass(i) * pl%vb(:, i) / cb%Gmass + if (pl%status(i) /= INACTIVE) cb%vb(:) = cb%vb(:) - pl%Gmass(i) * pl%vb(:, i) / cb%Gmass end do do concurrent(i = 1:npl) pl%vh(:, i) = pl%vb(:, i) - cb%vb(:) From e5b28e479aceb9f6f0c8149b58c05841e4dff954 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Mar 2023 11:19:23 -0400 Subject: [PATCH 7/7] Added the central body's contribution to orbital pre-collision angular momentum to the post-collision spin. --- src/symba/symba_discard.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index aa1f52296..a46ce057d 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -186,7 +186,7 @@ subroutine symba_discard_conserve_energy_and_momentum(pl, nbody_system, param, i becb1 = -(3 * cb%Gmass * cb%mass) / (5 * cb%radius) ! Add planet angular momentum to central body accumulator - cb%dL(:) = Lpl(:) + cb%dL(:) + cb%dL(:) = Lpl(:) + cb%dL(:) + Lcb(:) ! Update rotation of central body to by consistent with its angular momentum if (param%lrotation) then drot0(:) = cb%L0(:)/ (cb%Ip(3) * cb%mass * cb%radius**2)