Revert "Lots of updates to try to get DL/L0 under control"

This reverts commit 9918a0d.

src/main/swiftest_symba.f90

@@ -152,38 +152,35 @@ program swiftest_symba
       do while ((t < tstop) .and. ((ntp0 == 0) .or. (ntp > 0)))
          call symba_step_eucl(t, dt, param,npl,ntp,symba_plA, symba_tpA, nplplenc, npltpenc,&
                plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list)
-
+
+         ldiscard_pl = .false. 
+         ldiscard_tp = .false.
+         lfrag_add = .false.
          call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard_pl)
-         if (ldiscard_pl) then
-            if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot)
-            call symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
-            if (param%lenergy) then
-               call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot)
-               Ecollision = Eorbit_before - Eorbit_after    ! Energy change resulting in this collisional event Total running energy offset from collision in this step
-               symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision
-            end if
-         end if
-
-         lfrag_add = any(plplenc_list%status(1:nplplenc) == COLLISION)
-         if (lfrag_add) then
-            if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot)
-            call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, param)
-            if (param%lenergy) then
-               call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot)
-               Ecollision = Eorbit_before - Eorbit_after    ! Energy change resulting in this collisional event Total running energy offset from collision in this step
-               symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision
-               !write(*,*) 'Dke_orbit: ',ke_orbit_before - ke_orbit_after
-               !write(*,*) 'Dke_spin:  ',ke_spin_before - ke_spin_after
-               !write(*,*) 'Dpe:       ',pe_before - pe_after
-            end if
-         end if
-
          call symba_discard_tp(t, npl, ntp, symba_plA, symba_tpA, dt, rmin, rmax, rmaxu, qmin, qmin_coord, &    
                                 qmin_alo, qmin_ahi, param%lrhill_present, ldiscard_tp)
+         call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, mergeadd_list, nmergeadd, param)
+         ldiscard_pl = ldiscard_pl .or. lfrag_add
 
-         if (ldiscard_pl .or. ldiscard_tp .or. lfrag_add) then
+         ! These next two blocks should be streamlined/improved but right now we treat discards separately from collisions so it has to be this way
+         if (ldiscard_pl .or. ldiscard_tp) then
+            if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot)
+            call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, &
+                               discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param)
             call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, &
                                         discard_plA%helio%swiftest, discard_tpA%helio%swiftest, mergeadd_list, mergesub_list, discard_file, param%lbig_discard) 
+            nmergeadd = 0
+            nmergesub = 0
+            nsppl = 0
+            nsptp = 0
+            nplm = count(symba_plA%helio%swiftest%mass(1:npl) > mtiny)
+
+            if (param%lenergy)  then
+               call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot)
+               Ecollision = Eorbit_before - Eorbit_after    ! Energy change resulting in this collisional event Total running energy offset from collision in this step
+               symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision
+            end if
+            !if (ntp > 0) call util_dist_index_pltp(nplm, ntp, symba_plA, symba_tpA)
          end if
 
          iloop = iloop + 1

src/modules/module_interfaces.f90

@@ -796,21 +796,20 @@ END SUBROUTINE symba_chk
      END INTERFACE
 
      INTERFACE
-      subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, param)
-         use swiftest_globals
-         use swiftest_data_structures
-         use module_helio
-         use module_symba
-         implicit none
-         real(DP), intent(in)                      :: t
-         type(symba_pl)                            :: symba_plA
-         integer(I4B), intent(inout)               :: nplplenc, nmergeadd
-         type(symba_plplenc), intent(inout)        :: plplenc_list
-         logical, intent(inout)                    :: lfrag_add
-         type(symba_merger), intent(inout)         :: mergeadd_list
-         type(symba_pl), intent(inout)             :: discard_plA
-         type(user_input_parameters),intent(inout) :: param
-      end subroutine symba_collision
+          SUBROUTINE symba_collision(t, symba_plA, nplplenc, plplenc_list, ldiscard, mergeadd_list, nmergeadd, param)
+               USE swiftest_globals
+               USE swiftest_data_structures
+               USE module_helio
+               USE module_symba
+               IMPLICIT NONE
+               real(DP), intent(in)                      :: t
+               integer(I4B), intent(inout)               :: nplplenc, nmergeadd
+               type(symba_pl)                            :: symba_plA
+               type(symba_plplenc), intent(inout)        :: plplenc_list
+               type(symba_merger), intent(inout)         :: mergeadd_list
+               logical, intent(inout)                    :: ldiscard
+               type(user_input_parameters),intent(inout) :: param
+          END SUBROUTINE symba_collision
      END INTERFACE
 
      INTERFACE
@@ -1075,19 +1074,24 @@ END SUBROUTINE symba_peri
      END INTERFACE
 
      INTERFACE
-      subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
+      subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, &
+                              discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param)
          use swiftest_globals
          use swiftest_data_structures
          use module_symba
          implicit none
          real(DP),                    intent(in)    :: t
-         integer(I4B),                intent(inout) :: npl
+         integer(I4B),                intent(inout) :: npl, nplm, ntp, nsppl, nsptp
          integer(I4B),                intent(in)    :: nmergeadd
          type(symba_pl),              intent(inout) :: symba_plA
+         type(symba_tp),              intent(inout) :: symba_tpA
+         type(symba_tp),              intent(inout) :: discard_tpA
          type(symba_pl),              intent(inout) :: discard_plA
          type(symba_merger),          intent(inout) :: mergeadd_list 
+         logical(LGT),                intent(in)    :: ldiscard_pl, ldiscard_tp 
+         real(DP),                    intent(in)    :: mtiny
          type(user_input_parameters), intent(in)    :: param
-      end subroutine symba_rearray_pl
+      end subroutine symba_rearray
 
      END INTERFACE  
 

src/symba/symba_collision.f90

@@ -1,4 +1,4 @@
-subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, param)
+subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, mergeadd_list, nmergeadd, param)
    !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton
    !!
    !! Check for merger between planets in SyMBA. If the user has turned on the FRAGMENTATION feature, it will call the 
@@ -12,16 +12,15 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer
    use module_symba
    use module_interfaces, EXCEPT_THIS_ONE => symba_collision
    implicit none
-   ! Arguments
+
    real(DP), intent(in)                      :: t
-   type(symba_pl)                            :: symba_plA
    integer(I4B), intent(inout)               :: nplplenc, nmergeadd
+   type(symba_pl)                            :: symba_plA
    type(symba_plplenc), intent(inout)        :: plplenc_list
-   logical, intent(inout)                    :: lfrag_add
    type(symba_merger), intent(inout)         :: mergeadd_list
-   type(symba_pl), intent(inout)             :: discard_plA
+   logical, intent(inout)                    :: ldiscard
    type(user_input_parameters),intent(inout) :: param
-   ! internals
+
    integer(I4B), parameter                 :: NRES = 3   !! Number of collisional product results
    integer(I4B)                            :: i, j, k, index_enc, index_coll, jtarg, jproj
    real(DP), dimension(NRES)               :: mass_res
@@ -46,11 +45,19 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer
    logical, dimension(nplplenc)               :: lplpl_collision
    logical, dimension(:), allocatable         :: lplpl_unique_parent
    integer(I4B), dimension(:), pointer        :: plparent  
-   logical                                    :: ldiscard
+
+   ! TESTING
+   logical, save   :: lfirst = .true.
+   real(DP), save  :: Minitial
+   real(DP)        :: Msystem, Madd, Mdiscard
 
    ! First determine the collisional regime for each colliding pair
    associate(npl => symba_plA%helio%swiftest%nbody, xbpl => symba_plA%helio%swiftest%xb, statpl => symba_plA%helio%swiftest%status, idpl => symba_plA%helio%swiftest%id, &
              idx1 => plplenc_list%index1, idx2 => plplenc_list%index2, plparent => symba_plA%kin%parent)
+      if (lfirst) then
+         Minitial = sum(symba_plA%helio%swiftest%mass(1:npl))
+         lfirst = .false.
+      end if 
       lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION
       ldiscard = any(lplpl_collision)
       if (.not.ldiscard) return
@@ -304,8 +311,6 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer
          end do
          deallocate(family)
 
-         call symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
-
       end do
    end associate
 

src/symba/symba_frag_pos.f90

@@ -110,15 +110,7 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
       !write(*,        "(' ---------------------------------------------------------------------------')")
       !write(*,fmtlabel) ' Q_loss      |',-Qloss / abs(Etot_before)
       !write(*,        "(' ---------------------------------------------------------------------------')")
-!
-      ! Put any residual angular momentum into the spin before constraining energy
-      !L_residual(:) = Ltot_before(:) - Ltot_after(:)
-      !L_spin_frag(:) = L_residual(:) / nfrag
-      !do i = 1, nfrag
-      !   ke_spin_after = ke_spin_after - m_frag(i) * Ip_frag(3,i) * rad_frag(i)**2 * norm2(rot_frag(:,i))
-      !   rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2)
-      !   ke_spin_after = ke_spin_after + m_frag(i) * Ip_frag(3,i) * rad_frag(i)**2 * norm2(rot_frag(:,i))
-      !end do
+
 
       ! Set the "target" ke_after (the value of the orbital kinetic energy that the fragments ought to have)
       ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss
@@ -140,17 +132,24 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
       Lmag_after = norm2(Ltot_after(:))
 
       lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP) 
+      if (.not.lmerge) then
+         L_residual(:) = Ltot_before(:) - Ltot_after(:)
+         L_spin_frag(:) = L_residual(:) / nfrag
+         do i = 1, nfrag
+            rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2)
+         end do
+      end if
 
-      !write(*,        "(' ---------------------------------------------------------------------------')")
-      !write(*,        "('             |    T_orb    T_spin         T         pe      Etot      Ltot')")
-      !write(*,        "(' ---------------------------------------------------------------------------')")
-      !write(*,fmtlabel) ' change      |',(ke_after - ke_before) / abs(Etot_before), &
-      !                                 (ke_spin_after - ke_spin_before)/ abs(Etot_before), &
-      !                                 (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
-      !                                 (pe_after - pe_before) / abs(Etot_before), &
-      !                                 (Etot_after - Etot_before) / abs(Etot_before), &
-      !                                 (Lmag_after - Lmag_before) / Lmag_before
-      !write(*,        "(' ---------------------------------------------------------------------------')")
+!      write(*,        "(' ---------------------------------------------------------------------------')")
+!      write(*,        "('             |    T_orb    T_spin         T         pe      Etot      Ltot')")
+!      write(*,        "(' ---------------------------------------------------------------------------')")
+!      write(*,fmtlabel) ' change      |',(ke_after - ke_before) / abs(Etot_before), &
+!                                       (ke_spin_after - ke_spin_before)/ abs(Etot_before), &
+!                                       (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
+!                                       (pe_after - pe_before) / abs(Etot_before), &
+!                                       (Etot_after - Etot_before) / abs(Etot_before), &
+!                                       (Lmag_after - Lma_before) / Lmag_before
+!      write(*,        "(' ---------------------------------------------------------------------------')")
       !****************************************************************************************************************
 
    end associate
@@ -230,10 +229,8 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L
          x_frag(:, i) = r_frag_norm * (frag_vec(1) * x_col_unit(:) + frag_vec(2) * y_col_unit(:))
       end do
       call symba_frag_pos_com_adjust(xcom, m_frag, x_frag)
-      v_frag(:,:) = 0.0_DP
 
       call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag)
-      call symba_frag_pos_com_adjust(vcom, m_frag, v_frag)
 
       return
    end subroutine symba_frag_pos_initialize_fragments
@@ -275,12 +272,13 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius
          rmag = norm2(x_frag(:,i))
          v_r_unit(:) = x_frag(:,i) / rmag
          call util_crossproduct(h_unit(:), v_r_unit(:), v_t_unit(:))  ! make a unit vector in the tangential velocity direction wrt the angular momentum plane
-         v_frag(:,i) = v_frag(:,i) + L_orb_frag_mag / (m_frag(i) * rmag) * v_t_unit(:) ! Distribute the angular momentum equally amongst the fragments
+         v_frag(:,i) = L_orb_frag_mag / (m_frag(i) * rmag) * v_t_unit(:) ! Distribute the angular momentum equally amongst the fragments
       end do
 
       return 
    end subroutine symba_frag_pos_ang_mtm
 
+
    subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_frag, ke_target, lmerge)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
@@ -336,6 +334,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_fr
          lmerge = .true.
       end if
 
+
       ! Shift the fragments into the system barycenter frame
       do i = 1, nfrag
          v_frag(:,i) = v_corrected * mtot / m_frag(i) * v_r_unit(:, i) + v_t(:, i)
@@ -445,13 +444,11 @@ subroutine symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_orbit, ke_spi
          call move_alloc(ltmp, lexclude)
       end if
 
-      where (lexclude(1:npl_new))
-         symba_plwksp%helio%swiftest%status(1:npl_new) = INACTIVE
-      elsewhere
-         symba_plwksp%helio%swiftest%status(1:npl_new) = ACTIVE
+      where (lexclude(1:npl))
+         symba_plwksp%helio%swiftest%status(1:npl) = INACTIVE
       end where
 
-      nplm = count((symba_plwksp%helio%swiftest%mass(:) > param%mtiny) .and. .not.lexclude(:))
+      nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny)
       call util_dist_index_plpl(npl_new, nplm, symba_plwksp)
       call symba_energy_eucl(npl_new, symba_plwksp, param%j2rp2, param%j4rp4, ke_orbit, ke_spin, pe, te, Ltot)