Merge branch 'debug'

python/swiftest/swiftest/io.py

@@ -129,7 +129,6 @@ def str2bool(input_str):
         raise ValueError(f"{input_str} is not recognized as boolean")
 
 
-
 def real2float(realstr):
     """
     Converts a Fortran-generated ASCII string of a real value into a numpy float type. Handles cases where double precision
@@ -199,6 +198,7 @@ def read_swiftest_param(param_file_name, param, verbose=True):
         print(f"{param_file_name} not found.")
     return param
 
+
 def reorder_dims(ds):
 
     # Re-order dimension coordinates so that they are in the same order as the Fortran side
@@ -212,6 +212,8 @@ def reorder_dims(ds):
     idx = {index_name: idx[index_name] for index_name in dim_order}
     ds = ds.reindex(idx)
     return ds
+
+
 def read_swifter_param(param_file_name, verbose=True):
     """
     Reads in a Swifter param.in file and saves it as a dictionary

python/swiftest/swiftest/simulation_class.py

@@ -342,7 +342,7 @@ def __init__(self,read_param: bool = False, read_data: bool = False, simdir: os.
         # If the file doesn't exist, flag it for now so we know to create it
         param_file_found = False
         if read_param or read_data:
-            if self.read_param(read_init_cond = not read_data):
+            if self.read_param(read_init_cond = read_data):
                 # We will add the parameter file to the kwarg list. This will keep the set_parameter method from
                 # overriding everything with defaults when there are no arguments passed to Simulation()
                 kwargs['param_file'] = self.param_file

src/collision/collision_module.f90

@@ -151,6 +151,7 @@ module collision
       integer(I4B)                            :: status              !! Status flag to pass to the collision list once the collision has been resolved
       integer(I4B)                            :: collision_id        !! ID number of this collision event
       integer(I4B)                            :: maxid_collision = 0 !! The current maximum collision id number
+      real(DP)                                :: min_mfrag           !! Minimum fragment mass
 
       ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1
       real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor

src/collision/collision_regime.f90

@@ -213,10 +213,11 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
       Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3)
       bcrit = rad1 / (rad1 + rad2)
       Qloss = 0.0_DP
-      U_binding = (3 * GC * Mtot**2) / (5 * Rp) ! LS12 eq. 27
-      ke = 0.5_DP * Mtot * Vimp**2
-      pe = - GC * m1 * m2 / norm2(rh2 - rh1)
-      Qmerge = ke + pe + U_binding 
+      ! Specific binding energy
+      U_binding = (3 * GC * Mtot) / (5 * Rp) ! LS12 eq. 27
+      ke = 0.5_DP * Vimp**2
+      pe = - GC * m1 * m2 / (Mtot * norm2(rh2 - rh1))
+      Qmerge = ke + pe + U_binding ! The specific energy lost if this were a perfect merger
 
       if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then 
          regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime
@@ -246,7 +247,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
                Mlr = m1
                Mslr = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag)
                regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run
-               Qloss = (c_star + 1.0_DP) * U_binding ! Qr 
+               Qloss = (c_star - 1.0_DP) * U_binding 
             end if 
          else if (Vimp > Verosion .and. Vimp < Vsupercat) then
             if (m2 < 0.001_DP * m1) then 
@@ -257,28 +258,31 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
             else 
                Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag)  ! LS12 eq (37)
                regime = COLLRESOLVE_REGIME_DISRUPTION !disruption
-               Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_erosion
+               Qloss = (c_star - 1.0_DP) * U_binding 
             end if 
          else if (Vimp > Vsupercat) then 
             Mlr = max(Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(ETA), min_mfrag)   !LS12 eq (44)
             Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag)  !LS12 eq (37)
             regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic
-            Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat
+            Qloss = (c_star - 1.0_DP) * U_binding 
          else 
             call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Error no regime found in symba_regime")
          end if 
       end if 
 
+      Mresidual = Mtot - Mlr - Mslr
+      if (Mresidual < 0.0_DP) then ! prevents final masses from going negative
+         Mlr = Mlr + Mresidual
+      end if
+
       if (Mslr > Mlr) then ! The second-largest fragment is actually larger than the largest, so we will swap them
          Mtmp = Mlr
          Mlr = Mslr
          Mslr = Mtmp
       end if
 
-      Mresidual = Mtot - Mlr - Mslr
-      if (Mresidual < 0.0_DP) then ! prevents final masses from going negative
-         Mlr = Mlr + Mresidual
-      end if
+      Qloss = Qloss * Mtot ! Convert specific energy loss to total energy loss in the system
+
 
       return 
 

src/collision/collision_util.f90

@@ -921,6 +921,8 @@ module subroutine collision_util_set_natural_scale_factors(self)
          fragments%Gmass(:)  = fragments%Gmass(:)  / (collider%dscale**3/collider%tscale**2)
          fragments%radius(:) = fragments%radius(:) / collider%dscale
          impactors%Qloss     = impactors%Qloss     / collider%Escale
+
+         collider%min_mfrag = collider%min_mfrag / collider%mscale
       end associate
 
       return
@@ -986,6 +988,7 @@ module subroutine collision_util_set_original_scale_factors(self)
          collider%pe(:)        = collider%pe(:)        * collider%Escale
          collider%be(:)        = collider%be(:)        * collider%Escale
          collider%te(:)        = collider%te(:)        * collider%Escale
+         collider%min_mfrag    = collider%min_mfrag    * collider%mscale
 
          collider%mscale = 1.0_DP
          collider%dscale = 1.0_DP

src/fraggle/fraggle_generate.f90

@@ -468,18 +468,19 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       class(swiftest_parameters),   intent(inout) :: param        !! Current run configuration parameters 
       logical,                      intent(out)   :: lfailure     !! Did the velocity computation fail?
       ! Internals
-      integer(I4B) :: i, j, loop, try, istart, nfrag, nlast, nsteps
+      integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps
       logical :: lhitandrun, lsupercat
       real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual
-      real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric
+      real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric, dM
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
-      real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove
+      real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove, volume
       ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have
-      real(DP)                :: vmin_guess = 1.1_DP 
-      real(DP)                :: vmax_guess = 10.0_DP
-      real(DP)                :: delta_v, volume
-      integer(I4B), parameter :: MAXLOOP = 200
-      integer(I4B), parameter :: MAXTRY = 1000
+      real(DP)                :: vmin_guess = 1.01_DP 
+      real(DP)                :: vmax_guess = 8.0_DP
+      real(DP)                :: delta_v, GC
+      integer(I4B), parameter :: MAXLOOP = 100
+      integer(I4B), parameter :: MAXTRY = 20
+      real(DP), parameter     :: mass_reduction_ratio = 0.1_DP ! Ratio of difference between first and second fragment mass to remove from the largest fragment in case of a failure
       real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP
       class(collision_fraggle), allocatable :: collider_local
       character(len=STRMAX) :: message
@@ -489,9 +490,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
          lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) 
          lsupercat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) 
 
+         GC = impactors%Gmass(1) / impactors%mass(1)
+
          allocate(collider_local, source=collider)
          associate(fragments => collider_local%fragments)
 
+            volume(:) = 4.0_DP / 3.0_DP * PI * (fragments%radius(:))**3
+            fragments%density(:) = fragments%mass(:) / volume(:)
+
             ! The fragments will be divided into two "clouds" based on identified origin body. 
             ! These clouds will collectively travel like two impactors bouncing off of each other. 
             where(fragments%origin_body(:) == 1)
@@ -517,8 +523,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
             E_residual_best = huge(1.0_DP)
             lfailure = .false.
             dE_metric = huge(1.0_DP)
+            dE_best = huge(1.0_DP)
 
-            outer: do try = 1, nfrag
+            outer: do try = 1, MAXTRY
                ! Scale the magnitude of the velocity by the distance from the impact point
                ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct  
                do concurrent(i = 1:nfrag)
@@ -582,7 +589,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   dE = collider_local%te(2) - collider_local%te(1) 
                   E_residual = dE + impactors%Qloss
 
-                  if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (E_residual_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity
+                  if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity
                      E_residual_best = E_residual
                      dE_best = dE
 
@@ -620,24 +627,23 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call fragments%set_coordinate_system()
 
                end do
-               if (dE_best < 0.0_DP) exit outer
                ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters
-               if (fragments%nbody == 2) exit outer
-               ! Reduce the number of fragments by one
-               nlast = fragments%nbody
-               fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast)
-               fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast)
-               fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1)
-               fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast)
-               volume = 4.0_DP / 3.0_DP * PI * ((fragments%radius(1))**3 + (fragments%radius(nlast))**3)
-               fragments%density(1) = fragments%mass(1) / volume
-               fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD)
-               fragments%Ip(:,nlast) = 0.0_DP
-               fragments%mass(nlast) = 0.0_DP
-               fragments%Gmass(nlast) = 0.0_DP
-               fragments%radius(nlast) = 0.0_DP
-               fragments%status(nlast) = INACTIVE
-               fragments%nbody = nlast - 1
+               ! Reduce the fragment masses and add it to the largest remenant and try again
+               if (any(fragments%mass(2:nfrag) > collider%min_mfrag)) then
+                  do i = 2, nfrag
+                     if (fragments%mass(i) > collider%min_mfrag) then
+                        dM = min(mass_reduction_ratio * fragments%mass(i), fragments%mass(i) - collider%min_mfrag)
+                        fragments%mass(i) = fragments%mass(i) - dM
+                        fragments%mass(1) = fragments%mass(1) + dM
+                     end if
+                  end do
+               else
+                  exit outer
+               end if
+               fragments%Gmass(:) = GC * fragments%mass(:)
+
+               volume(:) = fragments%mass(:) / fragments%density(:)
+               fragments%radius(:) = (3._DP * volume(:) / (4._DP * PI))**(THIRD)
 
                call fragments%reset()
                call fraggle_generate_pos_vec(collider_local)

src/fraggle/fraggle_util.f90

@@ -25,7 +25,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       integer(I4B)              :: i, j, jproj, jtarg, nfrag, istart
       real(DP), dimension(2)    :: volume
       real(DP), dimension(NDIM) :: Ip_avg
-      real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul, G
+      real(DP) :: mfrag, mremaining, mtot, mcumul, G
       real(DP), dimension(:), allocatable :: mass
       real(DP), parameter :: BETA = 2.85_DP
       integer(I4B), parameter :: NFRAGMAX = 100  !! Maximum number of fragments that can be generated
@@ -37,7 +37,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       integer(I4B), dimension(:), allocatable :: ind
       logical :: flipper
 
-      associate(impactors => self%impactors)
+      associate(impactors => self%impactors, min_mfrag => self%min_mfrag)
          ! Get mass weighted mean of Ip and density
          volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3
          mtot = sum(impactors%mass(:))

src/swiftest/swiftest_util.f90

@@ -1869,21 +1869,21 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
                call plplenc_old%copy(nbody_system%plpl_encounter)
             end if
 
-            ! ! Re-build the encounter list
-            ! ! Be sure to get the level info if this is a SyMBA nbody_system
-            ! select type(nbody_system)
-            ! class is (symba_nbody_system)
-            ! select type(pl)
-            ! class is (symba_pl)
-            ! select type(tp)
-            ! class is (symba_tp)
-            !    lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) 
-            !    if (tp%nbody > 0) then
-            !       lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
-            !    end if
-            ! end select
-            ! end select
-            ! end select
+            ! Re-build the encounter list
+            ! Be sure to get the level info if this is a SyMBA nbody_system
+            select type(nbody_system)
+            class is (symba_nbody_system)
+            select type(pl)
+            class is (symba_pl)
+            select type(tp)
+            class is (symba_tp)
+               lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) 
+               if (tp%nbody > 0) then
+                  lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
+               end if
+            end select
+            end select
+            end select
 
             ! Re-index the encounter list as the index values may have changed
             if (nenc_old > 0) then