Merge branch 'debug'

src/fragmentation/fragmentation.f90

@@ -98,6 +98,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       L_frag_budget(:) = L_frag_orb(:) + L_frag_spin(:)
       f_spin = norm2(L_frag_spin(:)) / norm2(L_frag_budget(:))
 
+      call reset_fragments()
       call define_coordinate_system()
       call construct_temporary_system()
 
@@ -110,7 +111,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       ke_tot_deficit = 0.0_DP
       do while (try < MAXTRY)
          lfailure = .false.
-         call reset_fragments()
 
          call set_fragment_position_vectors()
 
@@ -145,6 +145,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
 
          if (.not.lfailure) exit
          call restructure_failed_fragments()
+         call reset_fragments()
          try = try + 1
       end do
       call restore_scale_factors()
@@ -326,8 +327,8 @@ subroutine define_coordinate_system()
             ! The cross product of the y- by z-axis will give us the x-axis
             x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:)
 
+            if (.not.any(x_frag(:,:) > 0.0_DP)) return
             rmag(:) = .mag. x_frag(:,:)
-            if (.not.any(rmag(:) > 0.0_DP)) return
 
             call random_number(L_sigma(:,:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector,
                                              ! otherwise we can get an ill-conditioned system

src/helio/helio_setup.f90

@@ -0,0 +1,22 @@
+submodule(helio_classes) s_helio_setup
+   use swiftest
+contains
+
+   module subroutine helio_setup_initialize_system(self, param)
+      !! author: David A. Minton
+      !!
+      !! Initialize a Helio nbody system from files, converting all heliocentric quantities to barycentric.
+      !!
+      implicit none
+      ! Arguments
+      class(helio_nbody_system),  intent(inout) :: self   !! Helio nbody system object
+      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+
+      call whm_setup_initialize_system(self, param)
+      call self%pl%h2b(self%cb)
+      call self%tp%h2b(self%cb)
+
+      return
+   end subroutine helio_setup_initialize_system
+
+end submodule s_helio_setup

src/io/io.f90

@@ -40,7 +40,8 @@ module subroutine io_conservation_report(self, param, lterminal)
                open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg)
             end if
          end if
-         call pl%h2b(cb)
+         call pl%vb2vh(cb)
+         call pl%xh2xb(cb)
          call system%get_energy_and_momentum(param) 
          ke_orbit_now = system%ke_orbit
          ke_spin_now = system%ke_spin
@@ -1234,13 +1235,18 @@ module subroutine io_write_discard(self, param)
       character(*), parameter :: NPLFMT    = '(I8)'
       character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))'
       class(swiftest_body), allocatable :: pltemp
-      character(len=STRMAX)   :: errmsg
+      character(len=STRMAX)   :: errmsg, out_stat
 
       if (param%discard_out == "") return
 
       associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody)
          if (nsp == 0) return
-         select case(param%out_stat)
+         if (lfirst) then
+            out_stat = param%out_stat
+         else
+            out_stat = 'APPEND'
+         end if
+         select case(out_stat)
          case('APPEND')
             open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg)
          case('NEW', 'REPLACE', 'UNKNOWN')

src/kick/kick.f90

@@ -14,7 +14,7 @@ module pure subroutine kick_getacch_int_pl(self)
       class(swiftest_pl), intent(inout) :: self
       ! Internals
       integer(I4B)                      :: k
-      real(DP)                          :: rji2, irij3, faci, facj
+      real(DP)                          :: rji2, irij3, faci, facj, rlim2
       real(DP)                          :: dx, dy, dz
 
       associate(pl => self, npl => self%nbody, nplpl => self%nplpl)
@@ -25,15 +25,18 @@ module pure subroutine kick_getacch_int_pl(self)
                   dy = pl%xh(2, j) - pl%xh(2, i)
                   dz = pl%xh(3, j) - pl%xh(3, i)
                   rji2 = dx**2 + dy**2 + dz**2
-                  irij3 = 1.0_DP / (rji2 * sqrt(rji2))
-                  faci = pl%Gmass(i) * irij3
-                  facj = pl%Gmass(j) * irij3
-                  pl%ah(1, i) = pl%ah(1, i) + facj * dx
-                  pl%ah(2, i) = pl%ah(2, i) + facj * dy
-                  pl%ah(3, i) = pl%ah(3, i) + facj * dz
-                  pl%ah(1, j) = pl%ah(1, j) - faci * dx
-                  pl%ah(2, j) = pl%ah(2, j) - faci * dy
-                  pl%ah(3, j) = pl%ah(3, j) - faci * dz
+                  rlim2 = (pl%radius(i) + pl%radius(j))**2
+                  if (rji2 > rlim2) then
+                     irij3 = 1.0_DP / (rji2 * sqrt(rji2))
+                     faci = pl%Gmass(i) * irij3
+                     facj = pl%Gmass(j) * irij3
+                     pl%ah(1, i) = pl%ah(1, i) + facj * dx
+                     pl%ah(2, i) = pl%ah(2, i) + facj * dy
+                     pl%ah(3, i) = pl%ah(3, i) + facj * dz
+                     pl%ah(1, j) = pl%ah(1, j) - faci * dx
+                     pl%ah(2, j) = pl%ah(2, j) - faci * dy
+                     pl%ah(3, j) = pl%ah(3, j) - faci * dz
+                  end if
                end if
             end associate
          end do

src/modules/helio_classes.f90

@@ -15,7 +15,8 @@ module helio_classes
    !********************************************************************************************************************************
    type, extends(whm_nbody_system) :: helio_nbody_system
    contains
-      procedure :: step => helio_step_system  !! Advance the Helio nbody system forward in time by one step
+      procedure :: step       => helio_step_system             !! Advance the Helio nbody system forward in time by one step
+      procedure :: initialize => helio_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates
    end type helio_nbody_system
 
    !********************************************************************************************************************************
@@ -35,8 +36,6 @@ module helio_classes
    !! Helio massive body particle class
    type, extends(swiftest_pl) :: helio_pl
    contains
-      procedure :: vh2vb       => helio_coord_vh2vb_pl     !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
-      procedure :: vb2vh       => helio_coord_vb2vh_pl     !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
       procedure :: drift       => helio_drift_pl           !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure :: lindrift    => helio_drift_linear_pl    !! Method for linear drift of massive bodies due to barycentric momentum of Sun
       procedure :: accel_gr    => helio_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction
@@ -53,8 +52,6 @@ module helio_classes
    !! Helio test particle class
    type, extends(swiftest_tp) :: helio_tp
    contains
-      procedure :: vh2vb       => helio_coord_vh2vb_tp     !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
-      procedure :: vb2vh       => helio_coord_vb2vh_tp     !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
       procedure :: lindrift    => helio_drift_linear_tp    !! Method for linear drift of massive bodies due to barycentric momentum of Sun
       procedure :: drift       => helio_drift_tp           !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure :: accel_gr    => helio_gr_kick_getacch_tp !! Acceleration term arising from the post-Newtonian correction
@@ -65,32 +62,6 @@ module helio_classes
    end type helio_tp
 
    interface
-      module subroutine helio_coord_vb2vh_pl(self, cb)
-         use swiftest_classes, only : swiftest_cb
-         implicit none
-         class(helio_pl),              intent(inout) :: self !! Helio massive body object
-         class(swiftest_cb),           intent(inout) :: cb   !! Swiftest central body object
-      end subroutine helio_coord_vb2vh_pl
-
-      module subroutine helio_coord_vb2vh_tp(self, vbcb)
-         implicit none
-         class(helio_tp),              intent(inout) :: self !! Helio massive body object
-         real(DP), dimension(:),       intent(in)    :: vbcb  !! Barycentric velocity of the central body
-      end subroutine helio_coord_vb2vh_tp
-
-      module subroutine helio_coord_vh2vb_pl(self, cb)
-         use swiftest_classes, only : swiftest_cb
-         implicit none
-         class(helio_pl),        intent(inout) :: self !! Helio massive body object
-         class(swiftest_cb),     intent(inout) :: cb   !! Swiftest central body object
-      end subroutine helio_coord_vh2vb_pl
-
-      module subroutine helio_coord_vh2vb_tp(self, vbcb)
-         implicit none
-         class(helio_tp),        intent(inout) :: self !! Helio massive body object
-         real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
-      end subroutine helio_coord_vh2vb_tp
-
       module subroutine helio_drift_body(self, system, param, dt)
          use swiftest_classes, only : swiftest_body, swiftest_nbody_system, swiftest_parameters
          implicit none
@@ -206,6 +177,13 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg)
          logical,                      intent(in)    :: lbeg   !! Logical flag indicating whether this is the beginning of the half step or not. 
       end subroutine helio_kick_vb_tp
 
+      module subroutine helio_setup_initialize_system(self, param)
+         use swiftest_classes, only : swiftest_parameters
+         implicit none
+         class(helio_nbody_system),  intent(inout) :: self   !! Helio nbody system object
+         class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+      end subroutine helio_setup_initialize_system
+
       module subroutine helio_step_pl(self, system, param, t, dt)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none