Fixed issue in new kick algorithm when there are semi-interacting bodies

examples/Basic_Simulation/initial_conditions.py

@@ -38,6 +38,7 @@
 # Initialize the simulation object as a variable. Arguments may be defined here or through the sim.run() method.
 #sim = swiftest.Simulation(fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)
 sim = swiftest.Simulation()
+sim.clean()
 
 # Add the modern planets and the Sun using the JPL Horizons Database.
 sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"])
@@ -73,9 +74,10 @@
 capm_tp     = default_rng().uniform(0.0, 360.0, ntp)
 
 sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
-# Display the run configuration parameters.
-sim.write_param()
-sim.get_parameter()
 
 # Run the simulation. Arguments may be defined here or thorugh the swiftest.Simulation() method.
-sim.run(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=10)
+#sim.run(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=10)
+sim.set_parameter(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=0)
+# Display the run configuration parameters.
+sim.get_parameter()
+sim.save()

examples/Fragmentation/Fragmentation_Movie.py

@@ -127,7 +127,7 @@
 nfrag_reduction = {"disruption_headon" : 1.0,
          "disruption_off_axis"         : 1.0,
          "supercatastrophic_headon"    : 1.0,
-         "supercatastrophic_off_axis"  : 1.0,
+         "supercatastrophic_off_axis"  : 10.0,
          "hitandrun_disrupt"           : 1.0,
          "hitandrun_pure"              : 1.0,
          "merge"                       : 1.0,
@@ -357,7 +357,7 @@ def vec_props(self, c):
 
         # Set fragmentation parameters
         minimum_fragment_gmass = 0.01 * body_Gmass[style][1] 
-        gmtiny = 0.10 * body_Gmass[style][1] 
+        gmtiny = 0.50 * body_Gmass[style][1] 
         sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style])
         sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0)
 

src/swiftest/swiftest_kick.f90

@@ -171,11 +171,12 @@ module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmas
       ! Internals
       integer(I4B) :: i, j
       real(DP)     :: rji2, rlim2, fac, rx, ry, rz
+      logical :: lmtiny
 
+      lmtiny = (nplm < npl)
       if (present(radius)) then
          !$omp parallel do default(private) schedule(static)&
-         !$omp shared(npl, nplm, r, Gmass, radius) &
-         !$omp reduction(-:acc) 
+         !$omp shared(npl,nplm, r, Gmass, radius, acc) 
          do i = 1, nplm
             do concurrent(j = 1:npl, j/=i)
                rx = r(1,j) - r(1,i)
@@ -184,31 +185,69 @@ module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmas
                rji2 = rx**2 + ry**2 + rz**2
                rlim2 = (radius(i) + radius(j))**2
                if (rji2 > rlim2)  then
-                  fac = Gmass(i) / (rji2 * sqrt(rji2))
-                  acc(1,j) = acc(1,j) - fac * rx
-                  acc(2,j) = acc(2,j) - fac * ry
-                  acc(3,j) = acc(3,j) - fac * rz
+                  fac = Gmass(j) / (rji2 * sqrt(rji2))
+                  acc(1,i) = acc(1,i) + fac * rx
+                  acc(2,i) = acc(2,i) + fac * ry
+                  acc(3,i) = acc(3,i) + fac * rz
                end if
             end do
          end do
          !$omp end parallel do
+
+         if (lmtiny) then
+            !$omp parallel do default(private) schedule(static)&
+            !$omp shared(npl,nplm, r, Gmass, radius, acc) 
+            do i = nplm+1,npl
+               do concurrent(j = 1:nplm, j/=i)
+                  rx = r(1,j) - r(1,i)
+                  ry = r(2,j) - r(2,i)
+                  rz = r(3,j) - r(3,i)
+                  rji2 = rx**2 + ry**2 + rz**2
+                  rlim2 = (radius(i) + radius(j))**2
+                  if (rji2 > rlim2)  then
+                     fac = Gmass(j) / (rji2 * sqrt(rji2))
+                     acc(1,i) = acc(1,i) + fac * rx
+                     acc(2,i) = acc(2,i) + fac * ry
+                     acc(3,i) = acc(3,i) + fac * rz
+                  end if
+               end do
+            end do
+            !$omp end parallel do
+         end if
       else
          !$omp parallel do default(private) schedule(static)&
-         !$omp shared(npl, nplm, r, Gmass) &
-         !$omp reduction(-:acc) 
+         !$omp shared(npl,nplm, r, Gmass, acc) 
          do i = 1, nplm
             do concurrent(j = 1:npl, j/=i)
                rx = r(1,j) - r(1,i)
                ry = r(2,j) - r(2,i)
                rz = r(3,j) - r(3,i)
-               rji2 = rx**2 + ry**2 + rz**2 
-               fac = Gmass(i) / (rji2 * sqrt(rji2))
-               acc(1,j) = acc(1,j) - fac * rx
-               acc(2,j) = acc(2,j) - fac * ry
-               acc(3,j) = acc(3,j) - fac * rz
+               rji2 = rx**2 + ry**2 + rz**2
+               fac = Gmass(j) / (rji2 * sqrt(rji2))
+               acc(1,i) = acc(1,i) + fac * rx
+               acc(2,i) = acc(2,i) + fac * ry
+               acc(3,i) = acc(3,i) + fac * rz
             end do
          end do
          !$omp end parallel do
+
+         if (lmtiny) then
+            !$omp parallel do default(private) schedule(static)&
+            !$omp shared(npl,nplm, r, Gmass, acc) 
+            do i = nplm+1,npl
+               do concurrent(j = 1:nplm, j/=i)
+                  rx = r(1,j) - r(1,i)
+                  ry = r(2,j) - r(2,i)
+                  rz = r(3,j) - r(3,i)
+                  rji2 = rx**2 + ry**2 + rz**2
+                  fac = Gmass(j) / (rji2 * sqrt(rji2))
+                  acc(1,i) = acc(1,i) + fac * rx
+                  acc(2,i) = acc(2,i) + fac * ry
+                  acc(3,i) = acc(3,i) + fac * rz
+               end do
+            end do
+            !$omp end parallel do
+         end if
       end if
 
       return