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/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() 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_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/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) 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 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(:) 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) diff --git a/src/walltime/walltime_implementations.f90 b/src/walltime/walltime_implementations.f90 index dcf1a2a60..0a9824d78 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,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 @@ -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 @@ -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 @@ -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 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)