From 15e766f078f7d740aaa00e469b852143c2ebdae9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 25 Feb 2023 16:22:30 -0500 Subject: [PATCH] Added more robust Brent's method solver borrowed and modified from CTEM --- src/fraggle/fraggle_util.f90 | 67 ++++-------- src/misc/solver_module.f90 | 191 ++++++++++++++++++++++++++++++++++- 2 files changed, 209 insertions(+), 49 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 3b7ee4fe2..14a3b69df 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -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) @@ -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 \ No newline at end of file diff --git a/src/misc/solver_module.f90 b/src/misc/solver_module.f90 index d0aa38e62..5bb3d1e0f 100644 --- a/src/misc/solver_module.f90 +++ b/src/misc/solver_module.f90 @@ -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 @@ -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) 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 \ No newline at end of file