Updated array bounds and initial conditions generator

MintonGroup · Jul 7, 2021 · ae446c9 · ae446c9
1 parent 886b63b
commit ae446c9
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 26 deletions.
diff --git a/examples/helio_swifter_comparison/init_cond.py b/examples/helio_swifter_comparison/init_cond.py
@@ -201,17 +201,17 @@
 
    sys.stdout = sys.__stdout__
    #Now make Swiftest files
-   #cbfile = open(swiftest_cb, 'w')
-   cbfile = FortranFile(swiftest_cb, 'w')
-   #print(1.0,file=cbfile)
-   #print(rmin,file=cbfile)
-   #print(J2,file=cbfile)
-   #print(J4,file=cbfile)
+   cbfile = open(swiftest_cb, 'w')
+   #cbfile = FortranFile(swiftest_cb, 'w')
+   print(1.0,file=cbfile)
+   print(rmin,file=cbfile)
+   print(J2,file=cbfile)
+   print(J4,file=cbfile)
    Msun = np.double(1.0)
-   cbfile.write_record(np.double(GMSun))
-   cbfile.write_record(np.double(rmin))
-   cbfile.write_record(np.double(J2))
-   cbfile.write_record(np.double(J4))
+   #cbfile.write_record(np.double(GMSun))
+   #cbfile.write_record(np.double(rmin))
+   #cbfile.write_record(np.double(J2))
+   #cbfile.write_record(np.double(J4))
    cbfile.close()
 
    #plfile = open(swiftest_pl, 'w')

diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90
@@ -69,25 +69,18 @@ module subroutine helio_drift_linear_pl(self, system, dt, pt)
       class(helio_nbody_system), intent(in)    :: system !! Swiftest nbody system object
       real(DP),                  intent(in)    :: dt     !! Stepsize
       real(DP), dimension(:),    intent(out)   :: pt     !! negative barycentric velocity of the central body
-      ! Internals
-      integer(I4B)          :: i
-      real(DP),dimension(NDIM) :: pttmp !intent(out) variables don't play nicely 
-                                        !with openmp's reduction for some reason
 
       associate(pl => self, npl => self%nbody, cb => system%cb)
-         pttmp(:) = 0.0_DP
-         do i = 2, npl
-            pttmp(:) = pttmp(:) + pl%Gmass(i) * pl%vb(:,i)
-         end do
-         pttmp(:) = pttmp(:) / cb%Gmass
-         do i = 2, npl
-            pl%xh(:,i) = pl%xh(:,i) + pttmp(:) * dt
-         end do
-         pt(:) = pttmp(:)
+         pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl))
+         pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl))
+         pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl))
+         pt(:) = pt(:) / cb%Gmass
+         pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt
+         pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt
+         pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt
       end associate
 
       return
-
    end subroutine helio_drift_linear_pl
 
    module subroutine helio_drift_tp(self, system, param, dt)

diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90
@@ -22,7 +22,7 @@ module subroutine helio_getacch_pl(self, system, param, t)
       real(DP), dimension(:, :), allocatable, save :: xh_loc, aobl
 
       associate(cb => system%cb, pl => self, npl => self%nbody)
-         pl%ahi(:,2:npl) = 0.0_DP
+         pl%ahi(:,:) = 0.0_DP
          call helio_getacch_int_pl(pl, t)
          pl%ah(:,:) = pl%ahi(:,:)
          if (param%loblatecb) call pl%obl_acc(cb)
@@ -84,7 +84,7 @@ subroutine helio_getacch_int_pl(pl, t)
          real(DP), dimension(NDIM) :: dx
 
          associate(npl => pl%nbody)
-            do i = 2, npl - 1
+            do i = 1, npl - 1
                do j = i + 1, npl
                   dx(:) = pl%xh(:,j) - pl%xh(:,i)
                   rji2 = dot_product(dx(:), dx(:))