Reverted orbital element utility functions to element-wise scalar ope…

…rations and got rid of .cross. , .mag. operators. These ran slow when coarrays were used (for some reason I don't entirly understand)

src/base/base_module.f90

@@ -30,6 +30,11 @@ module base
       real(DP)                                :: dt                   = -1.0_DP         !! Time step
       integer(I8B)                            :: iloop                = 0_I8B           !! Main loop counter
       integer(I8B)                            :: nloops               = 0_I8B           !! Total number of loops to execute
+      integer(I8B)                            :: istart               = 0_I8B           !! Starting index for loop counter
+      integer(I4B)                            :: iout                 = 0               !! Output cadence counter
+      integer(I4B)                            :: idump                = 0               !! Dump cadence counter
+      integer(I4B)                            :: nout                 = 0               !! Current output step
+      integer(I4B)                            :: istep                = 0               !! Current value of istep (used for time stretching)
       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

src/collision/collision_module.f90

@@ -19,7 +19,11 @@ module collision
    public
 
    character(len=*), parameter :: COLLISION_OUTFILE = 'collisions.nc'  !! Name of NetCDF output file for collision information
+#ifdef COARRAY
+   character(len=STRMAX) :: COLLISION_LOG_OUT !! Name of log file for collision diagnostic information (each co-image gets its own)
+#else
    character(len=*), parameter :: COLLISION_LOG_OUT = "collisions.log" !! Name of log file for collision diagnostic information
+#endif
 
    !>Symbolic names for collisional outcomes from collresolve_resolve:
    integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE              =  1

src/rmvs/rmvs_discard.f90

@@ -54,6 +54,9 @@ module subroutine rmvs_discard_tp(self, nbody_system, param)
          call swiftest_discard_tp(tp, nbody_system, param)
       end associate
 
+
+      return
+
    end subroutine rmvs_discard_tp
 
 end submodule s_rmvs_discard

src/swiftest/swiftest_discard.f90

@@ -131,7 +131,7 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param)
       ! Internals
       integer(I4B)        :: i
       real(DP)            :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2
-      character(len=STRMAX) :: idstr, timestr
+      character(len=STRMAX) :: idstr, timestr, message
 
       associate(ntp => tp%nbody, cb => nbody_system%cb, Gmtot => nbody_system%Gmtot)
          rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius)
@@ -144,8 +144,9 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param)
                   tp%status(i) = DISCARDED_RMAX
                   write(idstr, *) tp%id(i)
                   write(timestr, *) nbody_system%t
-                  write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
+                  write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
                               " too far from the central body at t = " // trim(adjustl(timestr))
+                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
                   tp%ldiscard(i) = .true.
                   tp%lmask(i) = .false.
                   call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), &
@@ -154,8 +155,9 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param)
                   tp%status(i) = DISCARDED_RMIN
                   write(idstr, *) tp%id(i)
                   write(timestr, *) nbody_system%t
-                  write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
+                  write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
                               " too close to the central body at t = " // trim(adjustl(timestr))
+                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
                   tp%ldiscard(i) = .true.
                   tp%lmask(i) = .false.
                   call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), &
@@ -168,8 +170,9 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param)
                      tp%status(i) = DISCARDED_RMAXU
                      write(idstr, *) tp%id(i)
                      write(timestr, *) nbody_system%t
-                     write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
+                     write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
                                  " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
+                     call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
                      tp%ldiscard(i) = .true.
                      tp%lmask(i) = .false.
                      call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), &
@@ -254,7 +257,7 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param)
       integer(I4B)              :: i, j, isp
       real(DP)                  :: r2min, radius
       real(DP), dimension(NDIM) :: dx, dv
-      character(len=STRMAX) :: idstri, idstrj, timestr
+      character(len=STRMAX) :: idstri, idstrj, timestr, message
 
       associate(ntp => tp%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, t => nbody_system%t, dt => param%dt)
          do i = 1, ntp
@@ -271,9 +274,10 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param)
                      write(idstri, *) tp%id(i)
                      write(idstrj, *) pl%id(j)
                      write(timestr, *) nbody_system%t
-                     write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstri)) // ")" &
+                     write(message, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstri)) // ")" &
                                                   // "  too close to massive body " // trim(adjustl(pl%info(j)%name)) // " ("  // trim(adjustl(idstrj)) // ")" &
                                                   // " at t = " // trim(adjustl(timestr))
+                     call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
                      tp%ldiscard(i) = .true.
                      call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), &
                                                discard_vh=tp%vh(:,i), discard_body_id=pl%id(j))

src/swiftest/swiftest_drift.f90

@@ -1,10 +1,10 @@
-!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh
+!! Coryright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh
 !! This file is part of Swiftest.
 !! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License 
 !! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
 !! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
 !! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-!! You should have received a copy of the GNU General Public License along with Swiftest. 
+!! You should have received a cory of the GNU General Public License along with Swiftest. 
 !! If not, see: https://www.gnu.org/licenses. 
 
 submodule (swiftest) s_swiftest_drift
@@ -104,7 +104,7 @@ module subroutine swiftest_drift_all(mu, x, v, n, param, dt, lmask, iflag)
    end subroutine swiftest_drift_all
 
 
-   pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) 
+   pure elemental module subroutine swiftest_drift_one(mu, rx, ry, rz, vx, vy, vz, dt, iflag) 
       !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
       !!
       !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient
@@ -114,18 +114,18 @@ pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz,
       implicit none
       ! Arguments
       real(DP), intent(in)      :: mu                     !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift
-      real(DP), intent(inout)   :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift
+      real(DP), intent(inout)   :: rx, ry, rz, vx, vy, vz !! Position and velocity of body to drift
       real(DP), intent(in)      :: dt                     !! Step size
       integer(I4B), intent(out) :: iflag                  !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       ! Internals
       integer(I4B) :: i
       real(DP)   :: dttmp
 
-      call swiftest_drift_dan(mu, px, py, pz, vx, vy, vz, dt, iflag)
+      call swiftest_drift_dan(mu, rx, ry, rz, vx, vy, vz, dt, iflag)
       if (iflag /= 0) then
          dttmp = 0.1_DP * dt
          do i = 1, 10
-            call swiftest_drift_dan(mu, px, py, pz, vx, vy, vz, dttmp, iflag)
+            call swiftest_drift_dan(mu, rx, ry, rz, vx, vy, vz, dttmp, iflag)
             if (iflag /= 0) exit
          end do
       end if
@@ -134,7 +134,7 @@ pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz,
    end subroutine swiftest_drift_one
 
 
-   pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag)
+   pure subroutine swiftest_drift_dan(mu, rx0, ry0, rz0, vx0, vy0, vz0, dt0, iflag)
       !! author: David A. Minton
       !!
       !! Perform Kepler drift, solving Kepler's equation in appropriate variables
@@ -144,23 +144,21 @@ pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag)
       implicit none
       ! Arguments
       real(DP),     intent(in)    :: mu            !! G * (m1 + m2), G = gravitational constant, m1 = mass of central body, m2 = mass of body to drift
-      real(DP),     intent(inout) :: px0, py0, pz0 !! position of body to drift
+      real(DP),     intent(inout) :: rx0, ry0, rz0 !! position of body to drift
       real(DP),     intent(inout) :: vx0, vy0, vz0 !! velocity of body to drift     
       real(DP),     intent(in)    :: dt0           !! time step
       integer(I4B), intent(out)   :: iflag         !! error status flag for Kepler drift (0 = OK, nonzero = NO CONVERGENCE)
       ! Internals
-      real(DP)        :: dt, f, g, fdot, gdot, c1, c2, c3, u, alpha, fp, r0
+      real(DP)        :: rx, ry, rz, vx, vy, vz, dt
+      real(DP)        :: f, g, fdot, gdot, c1, c2, c3, u, alpha, fp, r0
       real(DP)        :: v0s, a, asq, en, dm, ec, es, esq, xkep, fchk, s, c
-      real(DP), dimension(NDIM) :: x, v, x0, v0
 
       ! Executable code
       iflag = 0
       dt = dt0
-      x0 = [px0, py0, pz0]
-      v0 = [vx0, vy0, vz0]
-      r0 = sqrt(dot_product(x0(:), x0(:))) 
-      v0s = dot_product(v0(:), v0(:))
-      u = dot_product(x0(:),  v0(:))
+      r0 = sqrt(rx0*rx0 + ry0*ry0 + rz0*rz0)
+      v0s = vx0*vx0 + vy0*vy0 + vz0*vz0
+      u = rx0*vx0 + ry0*vy0 + rz0*vz0
       alpha = 2 * mu / r0 - v0s
       if (alpha > 0.0_DP) then
          a = mu / alpha
@@ -186,10 +184,19 @@ pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag)
             g = dt + (s - xkep) / en
             fdot = -(a / (r0 * fp)) * en * s
             gdot = (c - 1.0_DP) / fp + 1.0_DP
-            x(:) = x0(:) * f + v0(:) * g
-            v(:) = x0(:) * fdot + v0(:) * gdot
-            px0 = x(1); py0 = x(2); pz0 = x(3)
-            vx0 = v(1); vy0 = v(2); vz0 = v(3)
+            rx = rx0 * f + vx0 * g
+            ry = ry0 * f + vy0 * g
+            rz = rz0 * f + vz0 * g
+            vx = rx0 * fdot + vx0 * gdot
+            vy = ry0 * fdot + vy0 * gdot
+            vz = rz0 * fdot + vz0 * gdot
+
+            rx0 = rx
+            ry0 = ry
+            rz0 = rz
+            vx0 = vx
+            vy0 = vy
+            vz0 = vz
             iflag = 0
             return
          end if
@@ -201,10 +208,19 @@ pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag)
          g = dt - mu * c3
          fdot = -mu / (fp * r0) * c1
          gdot = 1.0_DP - mu / fp * c2
-         x(:) = x0(:) * f + v0(:) * g
-         v(:) = x0(:) * fdot + v0(:) * gdot
-         px0 = x(1); py0 = x(2); pz0 = x(3)
-         vx0 = v(1); vy0 = v(2); vz0 = v(3)
+         rx = rx0 * f + vx0 * g
+         ry = ry0 * f + vy0 * g
+         rz = rz0 * f + vz0 * g
+         vx = rx0 * fdot + vx0 * gdot
+         vy = ry0 * fdot + vy0 * gdot
+         vz = rz0 * fdot + vz0 * gdot
+
+         rx0 = rx
+         ry0 = ry
+         rz0 = rz
+         vx0 = vx
+         vy0 = vy
+         vz0 = vz
       end if
 
       return

src/swiftest/swiftest_driver.f90

@@ -25,11 +25,6 @@ program swiftest_driver
    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
    character(len=:), allocatable             :: display_style     !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD"
-   integer(I8B)                              :: istart            !! Starting index for loop counter
-   integer(I4B)                              :: iout              !! Output cadence counter
-   integer(I4B)                              :: idump             !! Dump cadence counter
-   integer(I4B)                              :: nout              !! Current output step
-   integer(I4B)                              :: istep             !! Current value of istep (used for time stretching)
    type(walltimer)                           :: integration_timer !! Object used for computing elapsed wall time
 
    call swiftest_io_get_args(integrator, param_file_name, display_style)
@@ -45,6 +40,11 @@ program swiftest_driver
       dt              => param%dt, &
       tstop           => param%tstop, &
       iloop           => param%iloop, &
+      istart          => param%istart, &
+      iout            => param%iout, &
+      idump           => param%idump, &
+      nout            => param%nout, &
+      istep           => param%istep, &
       nloops          => param%nloops, &
       istep_out       => param%istep_out, &
       fstep_out       => param%fstep_out, &

src/swiftest/swiftest_io.f90

@@ -2300,9 +2300,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
          ! Calculate the G for the nbody_system units
          param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2))
 
-         ! A minimal log of collision outcomes is stored in the following log file
-         ! More complete data on collisions is stored in the NetCDF output files
-         call swiftest_io_log_start(param, COLLISION_LOG_OUT, "Collision logfile")
 
          if ((param%encounter_save /= "NONE")       .and. &
              (param%encounter_save /= "TRAJECTORY") .and. &
@@ -2451,7 +2448,11 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
                if (param%log_output) flush(param%display_unit) 
 #ifdef COARRAY
             end if !(this_image() == 1)
+            write(COLLISION_LOG_OUT,'("collision_coimage",I0.3,".log")') this_image()
 #endif
+            ! A minimal log of collision outcomes is stored in the following log file
+            ! More complete data on collisions is stored in the NetCDF output files
+            call swiftest_io_log_start(param, COLLISION_LOG_OUT, "Collision logfile")
          end if
          ! Print the contents of the parameter file to standard output
       end select

src/swiftest/swiftest_module.f90

@@ -521,11 +521,11 @@ module subroutine swiftest_drift_body(self, nbody_system, param, dt)
          real(DP),                          intent(in)    :: dt     !! Stepsize
       end subroutine swiftest_drift_body
 
-      pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag)
+      pure elemental module subroutine swiftest_drift_one(mu, rx, ry, rz, vx, vy, vz, dt, iflag)
          !$omp declare simd(swiftest_drift_one)
          implicit none
          real(DP),     intent(in)       :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift
-         real(DP),     intent(inout)    :: px, py, pz, vx, vy, vz  !! Position and velocity of body to drift
+         real(DP),     intent(inout)    :: rx, ry, rz, vx, vy, vz  !! Position and velocity of body to drift
          real(DP),     intent(in)       :: dt    !! Step size
          integer(I4B), intent(out)      :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       end subroutine swiftest_drift_one
@@ -1027,33 +1027,33 @@ pure module subroutine swiftest_orbel_scget(angle, sx, cx)
          real(DP), intent(out) :: sx, cx
       end subroutine swiftest_orbel_scget
 
-      pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
+      pure elemental module subroutine swiftest_orbel_xv2aeq(mu, rx, ry, rz, vx, vy, vz, a, e, q)
          !$omp declare simd(swiftest_orbel_xv2aeq)
          implicit none
          real(DP), intent(in)  :: mu       !! Gravitational constant
-         real(DP), intent(in)  :: px,py,pz !! Position vector
+         real(DP), intent(in)  :: rx,ry,rz !! Position vector
          real(DP), intent(in)  :: vx,vy,vz !! Velocity vector
          real(DP), intent(out) :: a        !! semimajor axis
          real(DP), intent(out) :: e        !! eccentricity
          real(DP), intent(out) :: q        !! periapsis
       end subroutine swiftest_orbel_xv2aeq
 
-      pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi)
+      pure module subroutine swiftest_orbel_xv2aqt(mu, rx, ry, rz, vx, vy, vz, a, q, capm, tperi)
          !$omp declare simd(swiftest_orbel_xv2aqt)
          implicit none
          real(DP), intent(in)  :: mu       !! Gravitational constant
-         real(DP), intent(in)  :: px,py,pz !! Position vector
+         real(DP), intent(in)  :: rx,ry,rz !! Position vector
          real(DP), intent(in)  :: vx,vy,vz !! Velocity vector
          real(DP), intent(out) :: a        !! semimajor axis
          real(DP), intent(out) :: q        !! periapsis
          real(DP), intent(out) :: capm     !! mean anomaly
          real(DP), intent(out) :: tperi    !! time of pericenter passage
       end subroutine swiftest_orbel_xv2aqt
 
-      pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
+      pure module subroutine swiftest_orbel_xv2el(mu, rx, ry, rz, vx, vy, vz, a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
          implicit none
          real(DP), intent(in)  :: mu    !! Gravitational constant
-         real(DP), intent(in)  :: px,py,pz !! Position vector
+         real(DP), intent(in)  :: rx,ry,rz !! Position vector
          real(DP), intent(in)  :: vx,vy,vz !! Velocity vector
          real(DP), intent(out) :: a     !! semimajor axis
          real(DP), intent(out) :: e     !! eccentricity