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

Commit

Permalink
Added more robust Brent's method solver borrowed and modified from CTEM
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Feb 25, 2023
1 parent 75f936b commit 15e766f
Show file tree
Hide file tree
Showing 2 changed files with 209 additions and 49 deletions.
67 changes: 19 additions & 48 deletions src/fraggle/fraggle_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,54 +170,10 @@ module subroutine fraggle_util_set_mass_dist(self, param)
mremaining = impactors%mass_dist(iMrem)

! The mass will be distributed in a power law where N>M=(M/Mslr)**(-beta/3)
! Use Newton's method solver to get the logspace slope of the mass function
! Use Brent's method solver to get the logspace slope of the mass function
Mrat = (mremaining + Mslr) / Mslr
x0 = -3.0_DP
x1 = 3.0_DP
do j = 1, MAXLOOP
y0 = Mrat - 1.0_DP
y1 = Mrat - 1.0_DP
do i = iMrem, nfrag
y0 = y0 - (i-1)**(-x0/3.0_DP)
y1 = y1 - (i-1)**(-x1/3.0_DP)
end do
if (y0*y1 < 0.0_DP) exit
if (y0 > 0.0_DP) x0 = x0 * 1.6_DP
if (y1 < 0.0_DP) x1 = x1 * 1.6_DP
end do

! Find the mass scaling factor with bisection
do j = 1, MAXLOOP
beta = 0.5_DP * (x0 + x1)
ymid = Mrat - 1.0_DP
do i = iMrem, nfrag
ymid = ymid - (i-1)**(-beta/3.0_DP)
end do
if (abs((y0 - ymid)/y0) < 1e-12_DP) exit
if (y0 * ymid < 0.0_DP) then
x1 = beta
y1 = ymid
else if (y1 * ymid < 0.0_DP) then
x0 = beta
y0 = ymid
else
exit
end if
! Finish with Newton if we still haven't made it
if (j == MAXLOOP) then
do k = 1, MAXLOOP
y = Mrat - 1.0_DP
yp = 0.0_DP
do i = iMrem, nfrag
y = y - (i-1)**(-beta/3.0_DP)
if (i > 3) yp = yp - (i-1)**(-beta/3.0_DP) * log(real(i-1,kind=DP))
end do
eps = y/yp
if (abs(eps) < epsilon(1.0_DP) * beta) exit
beta = beta + eps
end do
end if
end do
beta = 3.0_DP
call solve_roots(sfd_function,beta)

do i = iMrem,nfrag
mass(i) = Mslr * (i-1)**(-beta/3.0_DP)
Expand Down Expand Up @@ -284,7 +240,22 @@ module subroutine fraggle_util_set_mass_dist(self, param)

call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)
return

contains

function sfd_function(x) result(y)
implicit none
real(DP), intent(in) :: x
real(DP) :: y
integer(I4B) :: i

y = Mrat - 1.0_DP
do i = iMrem, nfrag
y = y - (i-1)**(-x/3.0_DP)
end do
end function sfd_function

end subroutine fraggle_util_set_mass_dist


end submodule s_fraggle_util
end submodule s_fraggle_util
191 changes: 190 additions & 1 deletion src/misc/solver_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,17 @@ module solver
use lambda_function
use, intrinsic :: ieee_exceptions
private
public :: solve_linear_system, solve_rkf45
public :: solve_linear_system, solve_rkf45, solve_roots

interface solve_linear_system
module procedure solve_linear_system_dp
module procedure solve_linear_system_qp
end interface

interface solve_roots
module procedure solve_roots_dp
end interface

! interface
! module function solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
! implicit none
Expand Down Expand Up @@ -260,4 +264,189 @@ end function ge_wpp
! end function solve_rkf45


subroutine solve_roots_dp(f,x,tol,lerr)
!! author: David A. Minton
!!
!! Uses Brent's method to find the root of a function f
implicit none
! Arguments
interface
function f(x) result(y)
import DP
real(DP), intent(in) :: x
real(DP) :: y
end function f
end interface
real(DP), intent(inout) :: x !! Initial guess and also the final answer
real(DP), optional, intent(in) :: tol !! The relative tolerance on the solution
logical, optional, intent(out) :: lerr !! Returns .true. if a root was found, otherwise returns .false.
! Internals
real(DP),parameter :: TOL_DEFAULT = 1.e-7_DP
integer(I4B),parameter :: maxIterations=100
real(DP) :: valueAtRoot, x1, startx1, x2, factor, Tolerance
real(DP),parameter :: FPP = 1.e-11_DP
real(DP),parameter :: nearzero = 1.e-20_DP
real(DP) :: resultat,AA,BB,CC,DD,EE,FA,FB,FC,Tol1,PP,QQ,RR,SS,xm
integer(I4B) :: i,j, ev,br, niter, error
integer(I4B),parameter :: NTRY = 50
integer(I4B),parameter :: NBRACKET =20
real(DP),parameter :: FIRSTFACTOR = 1.6_DP
real(DP) :: f1, f2, fmin
integer(I4B) :: numev,numbr

if (present(tol)) then
Tolerance = tol
else
Tolerance = TOL_DEFAULT
end if

factor=FIRSTFACTOR
numev = 0
numbr = 0
startx1 = x
everything: do ev = 1, NTRY
numev = numev + 1
if (ev == NTRY) then
if (present(lerr)) lerr = .true.
return
end if
bracket: do br = 1, NBRACKET
numbr = numbr + 1
x1 = startx1
x2 = startx1 + abs(startx1 * factor)

! First bracket the root
f1 = f(x1)
f2 = f(x2)
fmin = abs(f1)
do j = 1, NTRY
if ((f1 * f2 < 0._DP)) exit bracket
if (abs(f1) < abs(f2)) then
x1 = x1 + factor * (x1 - x2)
f1 = f(x1)
if (abs(f1) < fmin) then
fmin = abs(f1)
startx1 = x1
end if
else
x2 = x2 + factor * (x2 - x1)
if (x2 < 0.0_DP) exit
f2 = f(x2)
if (abs(f2) < fmin) then
fmin = abs(f2)
startx1 = x2
end if
end if
end do
x1 = x2
factor = 0.5_DP * (factor + 1._DP)
end do bracket

! Now do a Brent's method to find the root
error = 0
AA = x1
BB = x2
FA = f1
FB = f2
CC = AA; FC = FA; DD = BB - AA; EE = DD
if (.not.RootBracketed(FA,FB)) then
error = -1
resultat = x1
else
FC = FB
do i = 1, maxIterations
if (.not.RootBracketed(FC,FB)) then
CC = AA; FC = FA; DD = BB - AA; EE = DD
end if
if (abs(FC) < abs(FB)) then
AA = BB; BB = CC; CC = AA
FA = FB; FB = FC; FC = FA
end if
Tol1 = 2 * FPP * abs(BB) + 0.5_DP * Tolerance
xm = 0.5_DP * (CC-BB)
if ((abs(xm) <= Tol1).or.(abs(FA) < nearzero)) then
! A root has been found
resultat = BB
valueAtRoot = f(resultat)
exit
else
if ((abs(EE) >= Tol1).and.(abs(FA) > abs(FB))) then
SS = FB / FA
if (abs(AA - CC) < nearzero) then
PP = 2 * xm * SS
QQ = 1._DP - SS
else
QQ = FA / FC
RR = FB / FC
PP = SS * (2 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1._DP))
QQ = (QQ - 1._DP) * (RR - 1._DP) * (SS - 1._DP)
end if
if (PP > nearzero) QQ = -QQ
PP = abs(PP)
if ((2*PP)<min(3*xm *QQ-abs(Tol1*QQ),abs(EE*QQ))) then
EE = DD
DD = PP / QQ
else
DD = xm
EE = DD
end if
else
DD = xm
EE = DD
end if
AA = BB
FA = FB
if (abs(DD) > Tol1) then
BB = BB + DD
else
if (xm > 0) then
BB = BB + abs(Tol1)
else
BB = BB - abs(Tol1)
end if
end if
FB = f(BB)
end if
end do
if (i >= maxIterations) error = -2
end if
niter = i
x = resultat
if (error == 0) exit
factor = 0.5_DP * (factor + 1.0_DP) ! Failed. Try again with a new factor
end do everything

if (present(lerr)) lerr = (error /= 0)

return

contains
! returns the minimum of two real numbers
real(DP) Function Minimum(x1,x2)
real(DP) x1,x2,resultat
if (x1 < x2) then
resultat = x1
else
resultat = x2
endif
Minimum = resultat
end function Minimum

! TRUE if x1*x2 negative
integer Function RootBracketed(x1,x2)
real(DP) x1,x2
integer resultat
if ((x1 > 0.and.x2 > 0).or.(x1 < 0.and.x2 < 0)) then
resultat = 0
else
resultat = 1
endif
RootBracketed = resultat
end function RootBracketed



end subroutine solve_roots_dp


end module solver

0 comments on commit 15e766f

Please sign in to comment.