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

src/main/swiftest_symba.f90

@@ -152,35 +152,38 @@ 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)
-         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
-
-         ! 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 (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(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_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
-            !if (ntp > 0) call util_dist_index_pltp(nplm, ntp, symba_plA, symba_tpA)
+         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)
+
+         if (ldiscard_pl .or. ldiscard_tp .or. lfrag_add) then
+            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) 
          end if
 
          iloop = iloop + 1

src/modules/module_interfaces.f90

@@ -796,20 +796,21 @@ END SUBROUTINE symba_chk
      END INTERFACE
 
      INTERFACE
-          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
+      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
      END INTERFACE
 
      INTERFACE
@@ -1074,24 +1075,19 @@ END SUBROUTINE symba_peri
      END INTERFACE
 
      INTERFACE
-      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)
+      subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
          use swiftest_globals
          use swiftest_data_structures
          use module_symba
          implicit none
          real(DP),                    intent(in)    :: t
-         integer(I4B),                intent(inout) :: npl, nplm, ntp, nsppl, nsptp
+         integer(I4B),                intent(inout) :: npl
          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
+      end subroutine symba_rearray_pl
 
      END INTERFACE  
 

src/symba/symba_collision.f90

@@ -1,4 +1,4 @@
-subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, mergeadd_list, nmergeadd, param)
+subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmergeadd, mergeadd_list, discard_plA, 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,15 +12,16 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg
    use module_symba
    use module_interfaces, EXCEPT_THIS_ONE => symba_collision
    implicit none
-
+   ! Arguments
    real(DP), intent(in)                      :: t
-   integer(I4B), intent(inout)               :: nplplenc, nmergeadd
    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
-   logical, intent(inout)                    :: ldiscard
+   type(symba_pl), intent(inout)             :: discard_plA
    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
@@ -45,19 +46,11 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg
    logical, dimension(nplplenc)               :: lplpl_collision
    logical, dimension(:), allocatable         :: lplpl_unique_parent
    integer(I4B), dimension(:), pointer        :: plparent  
-
-   ! TESTING
-   logical, save   :: lfirst = .true.
-   real(DP), save  :: Minitial
-   real(DP)        :: Msystem, Madd, Mdiscard
+   logical                                    :: ldiscard
 
    ! 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
@@ -311,6 +304,8 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg
          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,7 +110,15 @@ 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
@@ -132,24 +140,17 @@ 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 - Lma_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 - Lmag_before) / Lmag_before
+      !write(*,        "(' ---------------------------------------------------------------------------')")
       !****************************************************************************************************************
 
    end associate
@@ -229,8 +230,10 @@ 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
@@ -272,13 +275,12 @@ 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) = L_orb_frag_mag / (m_frag(i) * rmag) * v_t_unit(:) ! Distribute the angular momentum equally amongst the fragments
+         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
       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
       !!
@@ -334,7 +336,6 @@ 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)
@@ -444,11 +445,13 @@ 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))
-         symba_plwksp%helio%swiftest%status(1:npl) = INACTIVE
+      where (lexclude(1:npl_new))
+         symba_plwksp%helio%swiftest%status(1:npl_new) = INACTIVE
+      elsewhere
+         symba_plwksp%helio%swiftest%status(1:npl_new) = ACTIVE
       end where
 
-      nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny)
+      nplm = count((symba_plwksp%helio%swiftest%mass(:) > param%mtiny) .and. .not.lexclude(:))
       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)