Greatly improved robustness of the fragmentation code in symba_frag_pos. Reduced the number of failed initializations substantially

src/modules/lambda_function.f90

@@ -4,16 +4,29 @@ module lambda_function
    implicit none
 
    type, public :: lambda_obj 
+      !! Base class for an lambda function object. This object takes no additional arguments other than the dependent variable x, an array of real numbers
       procedure(lambda0), pointer, nopass :: lambdaptr => null()
-
+      real(DP) :: lastval
+      real(DP),dimension(:), allocatable :: lastarg
    contains
       generic   :: init => lambda_init_0
       procedure :: eval => lambda_eval_0
       procedure, nopass :: lambda_init_0
       final     :: lambda_destroy
    end type
+
+   type, public, extends(lambda_obj) :: lambda_obj_err
+      !! Extended class for an lambda function object. This object takes allows for the return of a logical error flag during evaluation of the function.
+      procedure(lambda0err), pointer, nopass :: lambdaptr_err => null()
+      logical   :: lerr     
+   contains
+      generic   :: init => lambda_init_0_err
+      procedure :: eval => lambda_eval_0_err
+      procedure, nopass :: lambda_init_0_err
+   end type
    interface lambda_obj
       module procedure lambda_init_0
+      module procedure lambda_init_0_err
    end interface
 
    abstract interface
@@ -23,6 +36,13 @@ function lambda0(x) result(y)
          real(DP), dimension(:), intent(in) :: x
          real(DP)                           :: y
       end function
+      function lambda0err(x, lerr) result(y)
+         ! Template for a 0 argument function that returns an error value
+         import DP
+         real(DP), dimension(:), intent(in)  :: x
+         logical,                intent(out) :: lerr
+         real(DP)                            :: y
+      end function
    end interface
 
    contains
@@ -34,22 +54,53 @@ type(lambda_obj) function lambda_init_0(lambda)
          lambda_init_0%lambdaptr => lambda
          return
       end function lambda_init_0
+
+      type(lambda_obj_err) function lambda_init_0_err(lambda, lerr)
+         implicit none
+         ! Arguments
+         procedure(lambda0err)  :: lambda
+         logical, intent(in) :: lerr
+         lambda_init_0_err%lambdaptr_err => lambda
+         lambda_init_0_err%lerr = lerr
+         return
+      end function lambda_init_0_err
 
       function lambda_eval_0(self, x) result(y)
          implicit none
          ! Arguments
-         class(lambda_obj),      intent(in) :: self
+         class(lambda_obj),      intent(inout) :: self
          real(DP), dimension(:), intent(in) :: x
          ! Result
          real(DP)                      :: y
 
          if (associated(self%lambdaptr)) then
             y = self%lambdaptr(x)
+            self%lastval = y
+            if (allocated(self%lastarg)) deallocate(self%lastarg)
+            allocate(self%lastarg, source=x)
          else
             error stop "Lambda function was not initialized"
          end if
       end function lambda_eval_0
 
+      function lambda_eval_0_err(self, x) result(y)
+         implicit none
+         ! Arguments
+         class(lambda_obj_err),  intent(inout) :: self
+         real(DP), dimension(:), intent(in) :: x
+         ! Result
+         real(DP)                      :: y
+
+         if (associated(self%lambdaptr_err)) then
+            y = self%lambdaptr_err(x, self%lerr)
+            self%lastval = y
+            if (allocated(self%lastarg)) deallocate(self%lastarg)
+            allocate(self%lastarg, source=x)
+         else
+            error stop "Lambda function was not initialized"
+         end if
+      end function lambda_eval_0_err
+
       subroutine lambda_destroy(self)
          implicit none
          type(lambda_obj) :: self

src/modules/module_interfaces.f90

@@ -1614,12 +1614,12 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
             use swiftest_globals
             use lambda_function
             implicit none
-            integer(I4B),           intent(in)  :: N
-            class(lambda_obj),      intent(in)  :: f
-            real(DP), dimension(:), intent(in)  :: x0_d
-            real(DP),               intent(in)  :: eps_d
-            logical,                intent(out) :: lerr
-            real(DP), dimension(:), allocatable :: x1_d
+            integer(I4B),           intent(in)    :: N
+            class(lambda_obj),      intent(inout) :: f
+            real(DP), dimension(:), intent(in)    :: x0_d
+            real(DP),               intent(in)    :: eps_d
+            logical,                intent(out)   :: lerr
+            real(DP), dimension(:), allocatable   :: x1_d
          end function util_minimize_bfgs
       end interface
 

src/symba/symba_collision.f90

@@ -34,7 +34,7 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg
    real(DP)                                :: mchild, mtot
    real(DP), dimension(NDIM)               :: xc, vc, xcom, vcom, xchild, vchild, xcrossv
    real(DP)                                :: mtiny_si
-   real(DP)                                :: mlr, mslr, msys, msys_new, Qloss
+   real(DP)                                :: mlr, mslr, msys, msys_new, Qloss, impact_parameter
    integer(I4B), dimension(:), allocatable :: family, collision_idx, unique_parent_idx
    integer(I4B)                            :: fam_size, istart, ncollisions, nunique_parent
    type family_array
@@ -153,6 +153,9 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg
             Ip(:, j) = mass(j) * symba_plA%helio%swiftest%Ip(:, idx_parent(j))
             ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip)
             L_spin(:, j)  = Ip(3, j) * radius(j)**2 * symba_plA%helio%swiftest%rot(:, idx_parent(j))
+
+            ! Use conservation of energy and angular momentum to adjust velocities and distances to be those equivalent to where they would be when the mutual radii are touching
+            !call vector_adjust(mass(j), x(:,j), v(:,j))
             if (nchild(j) > 0) then
                do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties
                   idx_child = parent_child_index_array(j)%idx(i + 1)
@@ -308,4 +311,33 @@ subroutine symba_collision (t, symba_plA, nplplenc, plplenc_list, ldiscard, merg
 
    return
 
+   contains
+      subroutine vector_adjust(mass, radius, x, v)
+         !! author: David A. Minton
+         !!
+         !! Uses conservation of angular momentum and energy to adjust the velocity vectors of two bodies to be what they would be if they were touching
+         implicit none
+         ! Arguments
+         real(DP), dimension(:), intent(in) :: mass, radius
+         real(DP), dimension(:), intent(inout) :: x,v
+         ! Internals
+         real(DP) :: mrat, vmag0, vmag, rmag0, rmag, B
+         real(DP), dimension(NDIM) :: h0
+
+         mrat = mass(1) * mass(2) / sum(mass(1:2))
+         vmag0 = norm2(v(:))
+         rmag0 = norm2(x(:))
+         rmag = sum(radius(1:2))
+
+         ! Solve vis-viva to get the new magnitude of the velocity
+         vmag = 2 * sqrt((0.5_DP * vmag0**2 * mrat + mass(1) * mass(2) * (1.0_DP / rmag - 1.0_DP / rmag0)) / mrat)
+
+         ! Use conservation of angular momentum to get the new impact parameter
+         call util_cross_product(x(:), v(:), h0(:))
+         B = norm2(h0) / (rmag * vmag)
+
+
+
+      end subroutine vector_adjust
+
 end subroutine symba_collision