Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Browse files Browse the repository at this point in the history
Fixed some minor issues to help with convergence and consistency
  • Loading branch information
daminton committed May 20, 2021
1 parent 0b4dd7d commit 4cfb29a
Show file tree
Hide file tree
Showing 2 changed files with 8 additions and 5 deletions.
5 changes: 2 additions & 3 deletions src/symba/symba_frag_pos.f90
Original file line number Diff line number Diff line change
Expand Up @@ -500,7 +500,7 @@ subroutine set_fragment_position_vectors()

! Place the fragments into a region that is big enough that we should usually not have overlapping bodies
! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down)
r_max = 1.5_DP
r_max = 1.0_DP

! We will treat the first fragment of the list as a special case. It gets initialized the maximum distance along the original impactor distance vector.
! This is done because in a regular disruption, the first body is the largest fragment.
Expand All @@ -522,9 +522,8 @@ subroutine set_fragment_position_vectors()
loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j)))
end do
end do
r_max = r_max + 0.2_DP
end do

x_frag(:, :) = x_frag(:, :) * r_max
call shift_vector_to_origin(m_frag, x_frag)

do i = 1, nfrag
Expand Down
8 changes: 6 additions & 2 deletions src/util/util_minimize_bfgs.f90
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
real(DP), dimension(:), allocatable :: x1_d
! 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
integer(I4B), parameter :: MAXLOOP = 500 !! Maximum number of loops before method is determined to have failed
real(QP), parameter :: gradeps = 1e-5_QP !! Tolerance for gradient calculations
real(QP), dimension(N) :: S !! Direction vectors
real(QP), dimension(N) :: Snorm !! normalized direction
Expand Down Expand Up @@ -110,13 +110,17 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d)
! update H matrix
H(:,:) = H(:,:) + ((1._QP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py
if (any(H(:,:) > sqrt(huge(1._QP)) / N)) then
lerr = .true.
write(*,*) 'BFGS did not converge after ',i,'iterations: H too big'
exit
end if
! 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) write(*,*) "BFGS ran out of loops!"
if (i == MAXLOOP) then
lerr = .true.
write(*,*) "BFGS ran out of loops!"
end if
end do
x1_d = x1
call ieee_get_flag(ieee_usual, fpe_flag)
Expand Down

0 comments on commit 4cfb29a

Please sign in to comment.