Removed intermediate diagnostics and re-wrote the BFGS minimizer to a…

…ccept maxloop as an argument
MintonGroup · Aug 16, 2021 · dc4a86a · dc4a86a
1 parent 4767b5e
commit dc4a86a
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 28 deletions.
diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90
@@ -127,18 +127,17 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
          if (.not.lfailure) call set_fragment_tan_vel(lfailure)
 
          if (lfailure) then
-            write(*,*) 'Failed to find tangential velocities'
-
+            !write(*,*) 'Failed to find tangential velocities'
          else 
             call set_fragment_radial_velocities(lfailure)
-            if (lfailure) write(*,*) 'Failed to find radial velocities'
+            !if (lfailure) write(*,*) 'Failed to find radial velocities'
             if (.not.lfailure) then
                call calculate_system_energy(linclude_fragments=.true.)
                if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then
-                  write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol
+                  !write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol
                   lfailure = .true.
                else if (abs(dLmag) / Lmag_before > Ltol) then
-                  write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before
+                  !write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before
                   lfailure = .true.
                end if
             end if
@@ -744,14 +743,14 @@ subroutine set_fragment_tan_vel(lerr)
 
             ! If we are over the energy budget, flag this as a failure so we can try again
             lerr = ((ke_frag_budget - ke_frag_spin - ke_frag_orbit) < 0.0_DP)
-            write(*,*) 'Tangential'
-            write(*,*) 'Failure? ',lerr
-            call calculate_fragment_ang_mtm()
-            write(*,*) '|L_remainder| : ',.mag.(L_frag_budget(:) - L_frag_tot(:)) / Lmag_before
-            write(*,*) 'ke_frag_budget: ',ke_frag_budget
-            write(*,*) 'ke_frag_spin  : ',ke_frag_spin
-            write(*,*) 'ke_tangential : ',ke_frag_orbit
-            write(*,*) 'ke_radial     : ',ke_frag_budget - ke_frag_spin - ke_frag_orbit
+            ! write(*,*) 'Tangential'
+            ! write(*,*) 'Failure? ',lerr
+            ! call calculate_fragment_ang_mtm()
+            ! write(*,*) '|L_remainder| : ',.mag.(L_frag_budget(:) - L_frag_tot(:)) / Lmag_before
+            ! write(*,*) 'ke_frag_budget: ',ke_frag_budget
+            ! write(*,*) 'ke_frag_spin  : ',ke_frag_spin
+            ! write(*,*) 'ke_tangential : ',ke_frag_orbit
+            ! write(*,*) 'ke_radial     : ',ke_frag_budget - ke_frag_spin - ke_frag_orbit
 
             return
          end subroutine set_fragment_tan_vel
@@ -882,12 +881,12 @@ subroutine set_fragment_radial_velocities(lerr)
             end do
             ke_frag_orbit = 0.5_DP * ke_frag_orbit
 
-            write(*,*) 'Radial'
-            write(*,*) 'Failure? ',lerr 
-            write(*,*) 'ke_frag_budget: ',ke_frag_budget
-            write(*,*) 'ke_frag_spin  : ',ke_frag_spin
-            write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
-            write(*,*) 'ke_remainder  : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin)
+            ! write(*,*) 'Radial'
+            ! write(*,*) 'Failure? ',lerr 
+            ! write(*,*) 'ke_frag_budget: ',ke_frag_budget
+            ! write(*,*) 'ke_frag_spin  : ',ke_frag_spin
+            ! write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
+            ! write(*,*) 'ke_remainder  : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin)
             lerr = .false.
 
             return

diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
@@ -1014,15 +1014,16 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
          real(DP),                     intent(in)    :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. 
       end subroutine util_rescale_system
 
-      module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1)
+      module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1)
          use lambda_function
          implicit none
          integer(I4B),           intent(in)    :: N
          class(lambda_obj),      intent(inout) :: f
          real(DP), dimension(:), intent(in)    :: x0
          real(DP),               intent(in)    :: eps
          logical,                intent(out)   :: lerr
-         real(DP), dimension(:), allocatable :: x1
+         integer(I4B),           intent(in)    :: maxloop
+         real(DP), dimension(:), allocatable   :: x1
       end function util_minimize_bfgs
 
       module subroutine util_peri_tp(self, system, param) 

diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90
@@ -1,15 +1,16 @@
 submodule (swiftest_classes) s_util_minimize_bfgs
    use swiftest
 contains  
-   module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1)
+   module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1)
       !! author: David A. Minton
       !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  - 
       !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables.  
       !! It recieves as input:
       !!   f%eval(x) : lambda function object containing the objective function as the eval metho
-      !!   N    :  Number of variables of function f
-      !!   x0   :  Initial starting value of x
-      !!   eps  :  Accuracy of 1 - dimensional minimization at each step
+      !!   N       : Number of variables of function f
+      !!   x0      : Initial starting value of x
+      !!   eps     : Accuracy of 1 - dimensional minimization at each step
+      !!   maxloop : Maximum number of loops to attempt to find a solution
       !! The outputs include
       !!   lerr :  Returns .true. if it could not find the minimum
       !! Returns
@@ -23,12 +24,12 @@ module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1)
       class(lambda_obj),      intent(inout) :: f
       real(DP), dimension(:), intent(in)    :: x0
       real(DP),               intent(in)    :: eps
+      integer(I4B),           intent(in)    :: maxloop
       logical,                intent(out)   :: lerr
       ! Result
       real(DP), dimension(:), allocatable :: x1
       ! Internals
       integer(I4B) ::  i, j, k, l, conv, num
-      integer(I4B), parameter :: MAXLOOP = 100 !! Maximum number of loops before method is determined to have failed 
       real(DP), parameter     :: graddelta = 1e-4_DP !! Delta x for gradient calculations
       real(DP), dimension(N) :: S               !! Direction vectors 
       real(DP), dimension(N,N) :: H             !! Approximated inverse Hessian matrix 
@@ -56,7 +57,7 @@ module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1)
          return
       end if
       grad1(:) = grad0(:)
-      do i = 1, MAXLOOP 
+      do i = 1, maxloop 
          !check for convergence
          conv = count(abs(grad1(:)) > eps)
          ! write(*,*) 'loop: ', i
@@ -114,7 +115,7 @@ module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1)
          ! Stop everything if there are any exceptions to allow the routine to fail gracefully
          call ieee_get_flag(ieee_usual, fpe_flag)
          if (any(fpe_flag)) exit 
-         if (i == MAXLOOP) then
+         if (i == maxloop) then
             lerr = .true.
             ! write(*,*) "BFGS ran out of loops!"
          end if