Merge remote-tracking branch 'origin/master'

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 

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()
 

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

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.)

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)

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

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

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(:)

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)  

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

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)