Merge branch 'debug'

Makefile

@@ -60,9 +60,9 @@ SWIFTEST_MODULES =   swiftest_globals.f90 \
 
 include Makefile.Defines
 
-MKL_ROOT = /apps/spack/bell/apps/intel-parallel-studio/cluster.2019.5-intel-19.0.5-4brgqlf/mkl/lib
-IMKL = -I$(MKLROOT)/include
-LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul 
+#MKL_ROOT = /apps/spack/bell/apps/intel-parallel-studio/cluster.2019.5-intel-19.0.5-4brgqlf/mkl/lib
+#IMKL = -I$(MKLROOT)/include
+#LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul 
 
 MODULES         = $(SWIFTEST_MODULES) $(USER_MODULES) 
 

Makefile.Defines

@@ -61,26 +61,27 @@ IPRODUCTION = -no-wrap-margin -O3 -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEA
 
 #gfortran flags
 GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none
-GPAR = -fopenmp #-ftree-parallelize-loops=4
+GPAR = -fopenmp -ftree-parallelize-loops=4
 GMEM =  -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
-GPRODUCTION = -O2 -ffree-line-length-none $(GPAR)
+GPRODUCTION = -O3 -ffree-line-length-none $(GPAR)
 
 #FFLAGS  = $(IDEBUG) $(SIMDVEC) $(PAR)
 #FFASTFLAGS  = $(IDEBUG) $(SIMDVEC) $(PAR)
 #FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT)  #$(ADVIXE_FLAGS)
 #FFLAGS = $(IPRODUCTION) -fp-model=fast $(OPTREPORT)  #$(ADVIXE_FLAGS)
 #FORTRAN = ifort
-AR     = xiar
+#AR     = xiar
 
 FORTRAN = gfortran
-#FFLAGS =  $(GDEBUG) # $(GMEM) $(GPAR)
-FFLAGS =  $(GPRODUCTION) #-fcheck=all #-Wall AR     = ar # DO NOT include in CFLAGS the "-c" option to compile object only
-FSTRICTFLAGS = $(FFLAGS)
+#FFLAGS =  $(GDEBUG)  $(GMEM) #$(GPAR)
+FFLAGS =  $(GPRODUCTION) #-g -fbacktrace -fcheck=all #-Wall 
+FSTRICTFLAGS=  $(GPRODUCTION) #-g -fbacktrace -fcheck=all #-Wall 
+AR     = ar # DO NOT include in CFLAGS the "-c" option to compile object only
 # this is done explicitly as needed in the Makefile
 
-CC              = icc
-#CC              = cc
+#CC              = icc
+CC              = cc
 CFLAGS          = -O3 -w -m64 -std=c99
 
 64_BIT_REALS    = -r8

src/drift/drift.f90

@@ -88,7 +88,6 @@ end subroutine drift_all
 
 
    module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) 
-      !$omp declare simd(drift_one)
       !! 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
@@ -119,7 +118,6 @@ end subroutine drift_one
 
 
    pure subroutine drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag)
-      !$omp declare simd(drift_dan)
       !! author: David A. Minton
       !!
       !! Perform Kepler drift, solving Kepler's equation in appropriate variables
@@ -197,7 +195,6 @@ end subroutine drift_dan
 
 
    pure subroutine drift_kepmd(dm, es, ec, x, s, c)
-      !$omp declare simd(drift_kepmd)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in difference form for an ellipse for small input dm and eccentricity
@@ -243,7 +240,6 @@ end subroutine drift_kepmd
 
 
    pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag)
-      !$omp declare simd(drift_kepu)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in universal variables
@@ -272,7 +268,6 @@ end subroutine drift_kepu
 
 
    pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f)
-      !$omp declare simd(drift_kepu_fchk)
       !! author: David A. Minton
       !!
       !! Computes the value of f, the function whose root we are trying to find in universal variables
@@ -303,7 +298,6 @@ end subroutine drift_kepu_fchk
 
 
    pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s)
-      !$omp declare simd(drift_kepu_guess)
       !! author: David A. Minton
       !!
       !! Compute initial guess for solving Kepler's equation using universal variables
@@ -348,7 +342,6 @@ end subroutine drift_kepu_guess
 
 
    pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag)
-      !$omp declare simd(drift_kepu_lag)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in universal variables using Laguerre's method
@@ -403,7 +396,6 @@ end subroutine drift_kepu_lag
 
 
    pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag)
-      !$omp declare simd(drift_kepu_new)
       !! author: David A. Minton
       !!
       !! Solve Kepler's equation in universal variables using Newton's method
@@ -455,7 +447,6 @@ end subroutine drift_kepu_new
 
 
    pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag)
-      !$omp declare simd(drift_kepu_p3solve)
       !! author: David A. Minton
       !!
       !! Computes real root of cubic involved in setting initial guess for solving Kepler's equation in universal variables
@@ -506,7 +497,6 @@ end subroutine drift_kepu_p3solve
 
 
    pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3)
-      !$omp declare simd(drift_kepu_stumpff)
       !! author: David A. Minton
       !!
       !! Compute Stumpff functions needed for Kepler drift in universal variables

src/encounter/encounter_check.f90

@@ -770,7 +770,6 @@ end subroutine encounter_check_all_triangular_pltp
 
 
    module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr)
-      !$omp declare simd(encounter_check_one)
       !! author: David A. Minton
       !!
       !! Determine whether a test particle and planet are having or will have an encounter within the next time step
@@ -838,6 +837,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in
       end associate
       if (nenc == 0) return
 
+
       allocate(index1(nenc))
       allocate(index2(nenc))
       if (present(lvdotr)) allocate(lvdotr(nenc))
@@ -868,17 +868,18 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in
    end subroutine encounter_check_collapse_ragged_list
 
 
-   pure subroutine encounter_check_make_ragged_list(lencounteri, ind_arr, lenci)
+   pure subroutine encounter_check_make_ragged_list(lencounteri, ind_arr, nenc,index2)
       implicit none
       ! Arguments
       logical, dimension(:), intent(in) :: lencounteri
       integer(I4B), dimension(:), intent(in) :: ind_arr
-      type(encounter_list), intent(out) :: lenci
+      integer(I4B),                            intent(out) :: nenc
+      integer(I4B), dimension(:), allocatable, intent(out) :: index2
 
-      lenci%nenc = count(lencounteri(:))
-      if (lenci%nenc > 0) then
-         allocate(lenci%index2(lenci%nenc))
-         lenci%index2(:) = pack(ind_arr(:), lencounteri(:)) 
+      nenc = count(lencounteri(:))
+      if (nenc > 0) then
+         allocate(index2(nenc))
+         index2(:) = pack(ind_arr(:), lencounteri(:)) 
       end if
 
       return
@@ -1017,7 +1018,7 @@ subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibegx, ie
             ibegyi = ibegy(i) 
             iendyi = iendy(i)
             call encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind(:), ibegxi, iendxi, ibegyi, iendyi, ibegx(:), iendx(:), ibegy(:), iendy(:), lencounteri(:))
-            call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i))
+            call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2)
          else
             lenc(i)%nenc = 0
          end if
@@ -1055,7 +1056,7 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibegx, iendx,
             ibegyi = ibegy(i) 
             iendyi = iendy(i)
             call encounter_check_sweep_aabb_one_single_list(n, ext_ind(:), ibegxi, iendxi, ibegyi, iendyi, ibegx(:), iendx(:), ibegy(:), iendy(:), lencounteri(:))
-            call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i))
+            call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2)
          else
             lenc(i)%nenc = 0
          end if
@@ -1067,7 +1068,6 @@ end subroutine encounter_check_sweep_aabb_all_single_list
 
 
    pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
-      !$omp declare simd(encounter_check_sweep_aabb_one_double_list)
       !! author: David A. Minton
       !!
       !! Performs a sweep operation on a single body. Encounters from the same lists not allowed (e.g. pl-tp encounters only)
@@ -1100,7 +1100,6 @@ end subroutine encounter_check_sweep_aabb_one_double_list
 
 
    pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
-      !$omp declare simd(encounter_check_sweep_aabb_one_single_list)
       !! author: David A. Minton
       !!
       !! Performs a sweep operation on a single body. Mutual encounters allowed (e.g. pl-pl)

src/fraggle/fraggle_generate.f90

@@ -547,21 +547,19 @@ function radial_objective_function(v_r_mag_input) result(fval)
             real(DP), dimension(frag%nbody)       :: kearr
             real(DP)                              :: keo, ke_radial, rotmag2, vmag2
 
-            associate(nfrag => frag%nbody)
-               allocate(v_shift, mold=frag%vb)
-               v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, frag%v_r_unit, frag%v_t_mag, frag%v_t_unit, frag%mass, frag%vbcom) 
-               !$omp do simd firstprivate(nfrag) lastprivate(rotmag2, vmag2) 
-               do i = 1,nfrag
-                  rotmag2 = frag%rot(1,i)**2 + frag%rot(2,i)**2 + frag%rot(3,i)**2
-                  vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2
-                  kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * rotmag2 + vmag2) 
-               end do
-               !$omp end do simd
-               keo = 2 * frag%ke_budget - sum(kearr(:))
-               ke_radial = frag%ke_budget - frag%ke_orbit - frag%ke_spin
-               ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for
-               fval = (keo / (2 * ke_radial))**2
-            end associate
+            allocate(v_shift, mold=frag%vb)
+            v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, frag%v_r_unit, frag%v_t_mag, frag%v_t_unit, frag%mass, frag%vbcom) 
+            !$omp do simd firstprivate(frag)
+            do i = 1,frag%nbody
+               rotmag2 = frag%rot(1,i)**2 + frag%rot(2,i)**2 + frag%rot(3,i)**2
+               vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2
+               kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * rotmag2 + vmag2) 
+            end do
+            !$omp end do simd
+            keo = 2 * frag%ke_budget - sum(kearr(:))
+            ke_radial = frag%ke_budget - frag%ke_orbit - frag%ke_spin
+            ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for
+            fval = (keo / (2 * ke_radial))**2
 
             return
          end function radial_objective_function

src/helio/helio_drift.f90

@@ -76,7 +76,6 @@ end subroutine helio_drift_tp
 
 
    pure elemental subroutine helio_drift_linear_one(xhx, xhy, xhz, ptx, pty, ptz, dt)
-      !$omp declare simd(helio_drift_linear_one)
       implicit none
       real(DP), intent(inout) :: xhx, xhy, xhz
       real(DP), intent(in) :: ptx, pty, ptz, dt

src/kick/kick.f90

@@ -108,7 +108,8 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad
          zr = x(3, j) - x(3, i) 
          rji2 = xr**2 + yr**2 + zr**2
          rlim2 = (radius(i) + radius(j))**2
-         if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
+         if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
+                                 ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
       end do
       !$omp end parallel do 
 
@@ -156,7 +157,8 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius
             zr = x(3, j) - x(3, i) 
             rji2 = xr**2 + yr**2 + zr**2
             rlim2 = (radius(i) + radius(j))**2
-            if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
+            if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
+                                    ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
          end do
       end do
       !$omp end parallel do
@@ -210,7 +212,6 @@ end subroutine kick_getacch_int_all_tp
 
 
    module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
-      !$omp declare simd(kick_getacch_int_one_pl)
       !! author: David A. Minton
       !!
       !! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies
@@ -242,7 +243,6 @@ end subroutine kick_getacch_int_one_pl
 
 
    module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az)
-      !$omp declare simd(kick_getacch_int_one_tp)
       !! author: David A. Minton
       !!
       !! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair.

src/modules/encounter_classes.f90

@@ -10,7 +10,7 @@ module encounter_classes
    integer(I4B), parameter :: SWEEPDIM = 2
 
    type :: encounter_list
-      integer(I4B)                              :: nenc   !! Total number of encounters
+      integer(I4B)                              :: nenc = 0  !! Total number of encounters
       logical,      dimension(:),   allocatable :: lvdotr !! relative vdotr flag
       integer(I4B), dimension(:),   allocatable :: status !! status of the interaction
       integer(I4B), dimension(:),   allocatable :: index1 !! position of the first body in the encounter
@@ -52,7 +52,7 @@ module encounter_classes
    interface
       module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, lencounter, lvdotr)
          implicit none
-         integer(I4B),                 intent(in)  :: nenc       !! Number of encounters in the encounter lists
+         integer(I4B),                 intent(in)  :: nenc         !! Number of encounters in the encounter lists
          integer(I4B), dimension(:),   intent(in)  :: index1     !! List of indices for body 1 in each encounter
          integer(I4B), dimension(:),   intent(in)  :: index2     !! List of indices for body 2 in each encounter1
          real(DP),     dimension(:,:), intent(in)  :: x1, v1     !! Array of indices of bodies 1
@@ -65,7 +65,7 @@ module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc
       end subroutine encounter_check_all
 
       module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, index1, index2, nenc)
-         import swiftest_parameters
+         use swiftest_classes, only: swiftest_parameters
          implicit none
          class(swiftest_parameters),              intent(inout) :: param  !! Current Swiftest run configuration parameter5s
          integer(I4B),                            intent(in)    :: npl    !! Total number of massive bodies
@@ -80,7 +80,7 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i
       end subroutine encounter_check_all_plpl
 
       module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc)
-         import swiftest_parameters
+         use swiftest_classes, only: swiftest_parameters
          implicit none
          class(swiftest_parameters),              intent(inout) :: param  !! Current Swiftest run configuration parameter5s
          integer(I4B),                            intent(in)    :: nplm   !! Total number of fully interacting massive bodies 
@@ -99,7 +99,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt,
       end subroutine encounter_check_all_plplm
 
       module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc)
-         import swiftest_parameters
+         use swiftest_classes, only: swiftest_parameters
          implicit none
          class(swiftest_parameters),              intent(inout) :: param  !! Current Swiftest run configuration parameter5s
          integer(I4B),                            intent(in)    :: npl    !! Total number of massive bodies 
@@ -175,6 +175,7 @@ module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radi
       end subroutine encounter_io_write_frame
 
       module subroutine encounter_io_write_list(self, pl, encbody, param)
+         use swiftest_classes, only : swiftest_pl, swiftest_body, swiftest_parameters
          implicit none
          class(encounter_list),  intent(in) :: self    !! Swiftest encounter list object
          class(swiftest_pl),         intent(in) :: pl      !! Swiftest massive body object

src/modules/rmvs_classes.f90

@@ -112,6 +112,7 @@ module subroutine rmvs_discard_tp(self, system, param)
       end subroutine rmvs_discard_tp
 
       module function rmvs_encounter_check_tp(self, param, system, dt) result(lencounter)
+         use swiftest_classes, only : swiftest_parameters
          implicit none
          class(rmvs_tp),             intent(inout) :: self       !! RMVS test particle object  
          class(swiftest_parameters), intent(inout) :: param      !! Current run configuration parameters 

src/operators/operator_cross.f90

@@ -8,7 +8,6 @@
 contains
 
    module pure function operator_cross_sp(A, B) result(C)
-      !$omp declare simd(operator_cross_sp)
       implicit none
       real(SP), dimension(:), intent(in) :: A, B
       real(SP), dimension(3) :: C
@@ -19,7 +18,6 @@ module pure function operator_cross_sp(A, B) result(C)
    end function operator_cross_sp
 
    module pure function operator_cross_dp(A, B) result(C)
-      !$omp declare simd(operator_cross_dp)
       implicit none
       real(DP), dimension(:), intent(in) :: A, B
       real(DP), dimension(3) :: C
@@ -30,7 +28,6 @@ module pure function operator_cross_dp(A, B) result(C)
    end function operator_cross_dp
 
    module pure function operator_cross_qp(A, B) result(C)
-      !$omp declare simd(operator_cross_qp)
       implicit none
       real(QP), dimension(:), intent(in) :: A, B
       real(QP), dimension(3) :: C
@@ -41,7 +38,6 @@ module pure function operator_cross_qp(A, B) result(C)
    end function operator_cross_qp
 
    module pure function operator_cross_i1b(A, B) result(C)
-      !$omp declare simd(operator_cross_i1b)
       implicit none
       integer(I1B), dimension(:), intent(in) :: A, B
       integer(I1B), dimension(3) :: C
@@ -52,7 +48,6 @@ module pure function operator_cross_i1b(A, B) result(C)
    end function operator_cross_i1b
 
    module pure function operator_cross_i2b(A, B) result(C)
-      !$omp declare simd(operator_cross_i2b)
       implicit none
       integer(I2B), dimension(:), intent(in) :: A, B
       integer(I2B), dimension(3) :: C
@@ -63,7 +58,6 @@ module pure function operator_cross_i2b(A, B) result(C)
    end function operator_cross_i2b
 
    module pure function operator_cross_i4b(A, B) result(C)
-      !$omp declare simd(operator_cross_i4b)
       implicit none
       integer(I4B), dimension(:), intent(in) :: A, B
       integer(I4B), dimension(3) :: C
@@ -74,7 +68,6 @@ module pure function operator_cross_i4b(A, B) result(C)
    end function operator_cross_i4b
 
    module pure function operator_cross_i8b(A, B) result(C)     
-      !$omp declare simd(operator_cross_i8b)
       implicit none
       integer(I8B), dimension(:), intent(in) :: A, B
       integer(I8B), dimension(3) :: C