diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 169c42e2c..b4d2eef1c 100644 --- 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 index ee2521916..686959d87 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1014,7 +1014,7 @@ 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 @@ -1022,7 +1022,8 @@ module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) 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 index 9a0e4e12e..01b57d868 100644 --- 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