Made numerous changes attempting to get momentum error under control.…

… Still not there yet

src/fragmentation/fragmentation.f90

@@ -49,7 +49,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       class(swiftest_nbody_system), allocatable :: tmpsys
       class(swiftest_parameters), allocatable   :: tmpparam
 
-
       if (nfrag < NFRAG_MIN) then
          write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given."
          lfailure = .true.
@@ -114,10 +113,9 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
             if (lfailure) write(*,*) 'Failed to find radial velocities'
             if (.not.lfailure) then
                call calculate_system_energy(linclude_fragments=.true.)
-
-               write(*,*) 'Qloss : ',Qloss
-               write(*,*) '-dEtot: ',-dEtot
-               write(*,*) 'delta : ',abs((dEtot + Qloss)) 
+               ! write(*,*) 'Qloss : ',Qloss
+               ! write(*,*) '-dEtot: ',-dEtot
+               ! write(*,*) 'delta : ',abs((dEtot + Qloss)) 
                if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then
                   write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol
                   lfailure = .true.
@@ -132,6 +130,10 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
          call restructure_failed_fragments()
          try = try + 1
       end do
+      call restore_scale_factors()
+      call calculate_system_energy(linclude_fragments=.true.)
+
+
       write(*,        "(' -------------------------------------------------------------------------------------')")
       write(*,        "('  Final diagnostic')")
       write(*,        "(' -------------------------------------------------------------------------------------')")
@@ -151,7 +153,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       end if
       write(*,        "(' -------------------------------------------------------------------------------------')")
 
-      call restore_scale_factors()
       call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Save the current halting modes so we can turn them off temporarily
 
       return 
@@ -248,6 +249,11 @@ subroutine restore_scale_factors()
             Ltot_after = Ltot_after * Lscale
             Lmag_after = Lmag_after * Lscale 
 
+            dLmag = norm2(Ltot_after(:) - Ltot_before(:)) 
+            dEtot = Etot_after - Etot_before
+
+            call tmpsys%rescale(tmpparam, mscale**(-1), dscale**(-1), tscale**(-1))
+
             mscale = 1.0_DP
             dscale = 1.0_DP
             vscale = 1.0_DP

src/modules/rmvs_classes.f90

@@ -163,15 +163,15 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask)
          implicit none
          class(rmvs_pl),                  intent(inout) :: self         !! RMVS massive body object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine rmvs_util_append_pl
 
       module subroutine rmvs_util_append_tp(self, source, lsource_mask)
          use swiftest_classes, only : swiftest_body
          implicit none
          class(rmvs_tp),                  intent(inout) :: self         !! RMVS test particle object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine rmvs_util_append_tp
 
       module subroutine rmvs_util_fill_pl(self, inserts, lfill_list)

src/modules/swiftest_classes.f90

@@ -836,35 +836,35 @@ module subroutine util_append_arr_char_string(arr, source, lsource_mask)
          implicit none
          character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr          !! Destination array 
          character(len=STRMAX), dimension(:), allocatable, intent(in)    :: source       !! Array to append 
-         logical,               dimension(:), optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical,               dimension(:),              intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_arr_char_string
 
       module subroutine util_append_arr_DP(arr, source, lsource_mask)
          implicit none
          real(DP), dimension(:), allocatable, intent(inout) :: arr          !! Destination array 
          real(DP), dimension(:), allocatable, intent(in)    :: source       !! Array to append 
-         logical,  dimension(:), optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical,  dimension(:),              intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_arr_DP
 
       module subroutine util_append_arr_DPvec(arr, source, lsource_mask)
          implicit none
          real(DP), dimension(:,:), allocatable, intent(inout) :: arr          !! Destination array 
          real(DP), dimension(:,:), allocatable, intent(in)    :: source       !! Array to append 
-         logical,  dimension(:),   optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical,  dimension(:),                intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_arr_DPvec
 
       module subroutine util_append_arr_I4B(arr, source, lsource_mask)
          implicit none
          integer(I4B), dimension(:), allocatable, intent(inout) :: arr          !! Destination array 
          integer(I4B), dimension(:), allocatable, intent(in)    :: source       !! Array to append 
-         logical,      dimension(:), optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical,      dimension(:),              intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_arr_I4B
 
       module subroutine util_append_arr_logical(arr, source, lsource_mask)
          implicit none
          logical, dimension(:), allocatable, intent(inout) :: arr          !! Destination array 
          logical, dimension(:), allocatable, intent(in)    :: source       !! Array to append 
-         logical, dimension(:), optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),              intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_arr_logical
    end interface
 
@@ -873,21 +873,21 @@ module subroutine util_append_body(self, source, lsource_mask)
          implicit none
          class(swiftest_body),            intent(inout) :: self   !! Swiftest body object
          class(swiftest_body),            intent(in)    :: source !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask  !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask  !! Logical mask indicating which elements to append to
       end subroutine util_append_body
 
       module subroutine util_append_pl(self, source, lsource_mask)
          implicit none
          class(swiftest_pl),              intent(inout) :: self         !! Swiftest massive body object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_pl
 
       module subroutine util_append_tp(self, source, lsource_mask)
          implicit none
          class(swiftest_tp),              intent(inout) :: self         !! Swiftest test particle object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine util_append_tp
 
       module subroutine util_coord_b2h_pl(self, cb)

src/modules/symba_classes.f90

@@ -492,14 +492,14 @@ module subroutine symba_util_append_arr_info(arr, source, lsource_mask)
          implicit none
          type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr          !! Destination array 
          type(symba_particle_info), dimension(:), allocatable, intent(in)    :: source       !! Array to append 
-         logical,                   dimension(:), optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical,                   dimension(:),              intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine symba_util_append_arr_info
 
       module subroutine symba_util_append_arr_kin(arr, source, lsource_mask)
          implicit none
          type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr          !! Destination array 
          type(symba_kinship), dimension(:), allocatable, intent(in)    :: source       !! Array to append 
-         logical,             dimension(:), optional,    intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical,             dimension(:),              intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine symba_util_append_arr_kin
    end interface
 
@@ -509,23 +509,23 @@ module subroutine symba_util_append_merger(self, source, lsource_mask)
          implicit none
          class(symba_merger),             intent(inout) :: self         !! SyMBA massive body object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine symba_util_append_merger
 
       module subroutine symba_util_append_pl(self, source, lsource_mask)
          use swiftest_classes, only : swiftest_body
          implicit none
          class(symba_pl),                 intent(inout) :: self         !! SyMBA massive body object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine symba_util_append_pl
 
       module subroutine symba_util_append_tp(self, source, lsource_mask)
          use swiftest_classes, only : swiftest_body
          implicit none
          class(symba_tp),                 intent(inout) :: self        !! SyMBA test particle object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine symba_util_append_tp
    end interface 
 

src/modules/whm_classes.f90

@@ -232,7 +232,7 @@ module subroutine whm_util_append_pl(self, source, lsource_mask)
          implicit none
          class(whm_pl),                   intent(inout) :: self         !! WHM massive body object
          class(swiftest_body),            intent(in)    :: source       !! Source object to append
-         logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+         logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
       end subroutine whm_util_append_pl
 
       module subroutine whm_util_spill_pl(self, discards, lspill_list, ldestructive)

src/rmvs/rmvs_util.f90

@@ -11,7 +11,7 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask)
       !! Arguments
       class(rmvs_pl),                  intent(inout) :: self         !! RMVS massive body object
       class(swiftest_body),            intent(in)    :: source       !! Source object to append
-      logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+      logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
 
       select type(source)
       class is (rmvs_pl)
@@ -44,7 +44,7 @@ module subroutine rmvs_util_append_tp(self, source, lsource_mask)
       !! Arguments
       class(rmvs_tp),                  intent(inout) :: self         !! RMVS test particle object
       class(swiftest_body),            intent(in)    :: source       !! Source object to append
-      logical, dimension(:), optional, intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+      logical, dimension(:),           intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
 
       select type(source)
       class is (rmvs_tp)

src/symba/symba_fragmentation.f90

@@ -133,11 +133,13 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
                      plnew%k2 = pl%k2(ibiggest)
                      plnew%tlag = pl%tlag(ibiggest)
                   end if
+                  call plnew%set_mu(cb)
+                  pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
 
                   ! Append the new merged body to the list and record how many we made
                   nstart = pl_adds%nbody + 1
                   nend = pl_adds%nbody + plnew%nbody
-                  call pl_adds%append(plnew)
+                  call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)])
                   pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                   call plnew%setup(0, param)
@@ -291,11 +293,13 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m
                      plnew%k2 = pl%k2(ibiggest)
                      plnew%tlag = pl%tlag(ibiggest)
                   end if
+                  call plnew%set_mu(cb)
+                  pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
 
                   ! Append the new merged body to the list and record how many we made
                   nstart = pl_adds%nbody + 1
                   nend = pl_adds%nbody + plnew%nbody
-                  call pl_adds%append(plnew)
+                  call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)])
                   pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                   call plnew%setup(0, param)
@@ -427,11 +431,13 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass,
                   plnew%k2 = pl%k2(ibiggest)
                   plnew%tlag = pl%tlag(ibiggest)
                end if
+               call plnew%set_mu(cb)
+               pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
 
                ! Append the new merged body to the list and record how many we made
                nstart = pl_adds%nbody + 1
                nend = pl_adds%nbody + plnew%nbody
-               call pl_adds%append(plnew)
+               call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)])
                pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                call plnew%setup(0, param)
@@ -570,11 +576,13 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family,
                      plnew%k2 = pl%k2(ibiggest)
                      plnew%tlag = pl%tlag(ibiggest)
                   end if
+                  call plnew%set_mu(cb)
+                  pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
 
                   ! Append the new merged body to the list and record how many we made
                   nstart = pl_adds%nbody + 1
                   nend = pl_adds%nbody + plnew%nbody
-                  call pl_adds%append(plnew)
+                  call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)])
                   pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                   call plnew%setup(0, param)