From 15e766f078f7d740aaa00e469b852143c2ebdae9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 25 Feb 2023 16:22:30 -0500 Subject: [PATCH 1/3] 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 From 65eaabf8d698fcb3622f4e23c35db0f3fbf241cb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 25 Feb 2023 17:39:26 -0500 Subject: [PATCH 2/3] Added ability to stretch time during output steps --- python/swiftest/swiftest/simulation_class.py | 34 +++++++++++++++++--- src/base/base_module.f90 | 2 ++ src/misc/solver_module.f90 | 6 ++-- src/swiftest/swiftest_io.f90 | 33 +++++++++++++++++++ 4 files changed, 67 insertions(+), 8 deletions(-) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 7285eeb51..c908d322d 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -111,14 +111,20 @@ def __init__(self,read_param: bool = False, istep_out : int, optional The number of time steps between output saves to file. *Note*: only `istep_out` or `tstep_out` can be set. Parameter input file equivalent: `ISTEP_OUT` + tstep_out : float, optional + The approximate time between when outputs are written to file. Passing this computes + `istep_out = floor(tstep_out/dt)`. *Note*: only `istep_out` or `tstep_out` can be set. + Parameter input file equivalent: None + nstep_out : int, optional + The total number of times that outputs are written to file. Passing this allows for a geometric progression of output steps: + `TSTART, f**0 * TSTEP_OUT, f**1 * TSTEP_OUT, f**2 * TSTEP_OUT, ..., f**(nstep_out-1) * TSTEP_OUT`, + where `f` is a factor that can stretch (or shrink) the time between outputs. Setting `nstep_out = int((tstart - tstop) / (tstep_out))` is + equivalent to the standard linear output (i.e. `f==1`) and is the same as not passing anything for this argument. + *Note*: Passing `nstep_out` requires passing either `istep_out` or `tstep_out` as well. dump_cadence : int, optional The number of output steps (given by `istep_out`) between when the saved data is dumped to a file. Setting it to 0 is equivalent to only dumping data to file at the end of the simulation. Default value is 10. Parameter input file equivalent: `DUMP_CADENCE` - tstep_out : float, optional - The approximate time between when outputs are written to file. Passing this computes - `istep_out = floor(tstep_out/dt)`. *Note*: only `istep_out` or `tstep_out` can be set. - Parameter input file equivalent: None init_cond_file_type : {"NETCDF_DOUBLE", "NETCDF_FLOAT", "ASCII"}, default "NETCDF_DOUBLE" The file type containing initial conditions for the simulation: * NETCDF_DOUBLE: A single initial conditions input file in NetCDF file format of type NETCDF_DOUBLE. @@ -571,6 +577,7 @@ def set_simulation_time(self, dt: float | None = None, istep_out: int | None = None, tstep_out: float | None = None, + nstep_out: int | None = None, dump_cadence: int | None = None, verbose: bool | None = None, **kwargs: Any @@ -588,10 +595,18 @@ def set_simulation_time(self, dt : float, optional The step size of the simulation. `dt` must be less than or equal to `tstop-dstart`. istep_out : int, optional - The number of time steps between outputs to file. *Note*: only `istep_out` or `tstep_out` can be set. + The number of time steps between output saves to file. *Note*: only `istep_out` or `tstep_out` can be set. + Parameter input file equivalent: `ISTEP_OUT` tstep_out : float, optional The approximate time between when outputs are written to file. Passing this computes `istep_out = floor(tstep_out/dt)`. *Note*: only `istep_out` or `tstep_out` can be set. + Parameter input file equivalent: None + nstep_out : int, optional + The total number of times that outputs are written to file. Passing this allows for a geometric progression of output steps: + `TSTART, f**0 * TSTEP_OUT, f**1 * TSTEP_OUT, f**2 * TSTEP_OUT, ..., f**(nstep_out-1) * TSTEP_OUT`, + where `f` is a factor that can stretch (or shrink) the time between outputs. Setting `nstep_out = int((tstart - tstop) / (tstep_out))` is + equivalent to the standard linear output (i.e. `f==1`) and is the same as not passing anything for this argument. + *Note*: Passing `nstep_out` requires passing either `istep_out` or `tstep_out` as well. dump_cadence : int, optional The number of output steps (given by `istep_out`) between when the saved data is dumped to a file. Setting it to 0 is equivalent to only dumping data to file at the end of the simulation. Default value is 10. @@ -672,6 +687,12 @@ def set_simulation_time(self, if istep_out is not None: self.param['ISTEP_OUT'] = int(istep_out) + + if nstep_out is not None: + if istep_out is None: + warnings.warn("nstep_out requires either istep_out or tstep_out to also be set", stacklevel=2) + else: + self.param['NSTEP_OUT'] = int(nstep_out) if dump_cadence is None: dump_cadence = self.param.pop("DUMP_CADENCE", 1) @@ -714,6 +735,7 @@ def get_simulation_time(self, arg_list: str | List[str] | None = None, verbose: "tstop": "TSTOP", "dt": "DT", "istep_out": "ISTEP_OUT", + "nstep_out": "NSTEP_OUT", "dump_cadence": "DUMP_CADENCE", } @@ -723,6 +745,7 @@ def get_simulation_time(self, arg_list: str | List[str] | None = None, verbose: "dt": self.TU_name, "tstep_out": self.TU_name, "istep_out": "", + "nstep_out": "", "dump_cadence": ""} tstep_out = None @@ -772,6 +795,7 @@ def set_parameter(self, verbose: bool = True, **kwargs): "dt": None, "istep_out": 1, "tstep_out": None, + "nstep_out": None, "dump_cadence": 10, "init_cond_file_type": "NETCDF_DOUBLE", "init_cond_file_name": None, diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 51d632625..6060d53a1 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -35,6 +35,8 @@ module base character(STRMAX) :: in_type = "NETCDF_DOUBLE" !! Data representation type of input data files character(STRMAX) :: in_form = "XV" !! Format of input data files ("EL" or ["XV"]) integer(I4B) :: istep_out = -1 !! Number of time steps between saved outputs + integer(I4B) :: nstep_out = -1 !! Total number of saved outputs + real(DP) :: fstep_out = 1.0_DP !! The output step time stretching factor character(STRMAX) :: outfile = BIN_OUTFILE !! Name of output binary file character(STRMAX) :: out_type = "NETCDF_DOUBLE" !! Binary format of output file character(STRMAX) :: out_form = "XVEL" !! Data to write to output file diff --git a/src/misc/solver_module.f90 b/src/misc/solver_module.f90 index 5bb3d1e0f..410116b2a 100644 --- a/src/misc/solver_module.f90 +++ b/src/misc/solver_module.f90 @@ -290,7 +290,7 @@ end function f 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),parameter :: FIRSTFACTOR = 1.1_DP real(DP) :: f1, f2, fmin integer(I4B) :: numev,numbr @@ -313,7 +313,7 @@ end function f bracket: do br = 1, NBRACKET numbr = numbr + 1 x1 = startx1 - x2 = startx1 + abs(startx1 * factor) + x2 = startx1 + abs(startx1 * (1.0_DP - factor)) ! First bracket the root f1 = f(x1) @@ -339,7 +339,7 @@ end function f end if end do x1 = x2 - factor = 0.5_DP * (factor + 1._DP) + factor = factor + 0.5_DP * (factor - 1._DP) end do bracket ! Now do a Brent's method to find the root diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index faf17424f..b1af394ec 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -1883,6 +1883,7 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i integer(I4B) :: nseeds, nseeds_from_file logical :: seed_set = .false. !! Is the random seed set in the input file? character(len=:), allocatable :: integrator + real(DP) :: tratio, y ! Parse the file line by line, extracting tokens then matching them up with known parameters if possible @@ -1931,6 +1932,8 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i param%in_form = param_value case ("ISTEP_OUT") read(param_value, *) param%istep_out + case ("NSTEP_OUT") + read(param_value, *) param%nstep_out case ("BIN_OUT") param%outfile = param_value case ("OUT_TYPE") @@ -2126,6 +2129,15 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i iostat = -1 return end if + if (param%nstep_out <= 0) then + param%nstep_out = int((param%tstop - param%t0) / (param%istep_out * param%dt)) + param%fstep_out = 1.0_DP ! Linear output time + else + param%fstep_out = 1._DP + tratio = (param%TSTOP - param%T0) / (param%istep_out * param%dt) + y = time_stretcher(param%fstep_out) + call solve_roots(time_stretcher,param%fstep_out) + end if if (param%dump_cadence < 0) then write(iomsg,*) 'Invalid DUMP_CADENCE. Must be a positive integer or 0.' iostat = -1 @@ -2332,6 +2344,27 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i return 667 continue write(*,*) "Error reading param file: ", trim(adjustl(iomsg)) + + contains + function time_stretcher(fstep_out) result(ans) + !! author: David A. Minton + !! + !! Equation for the time stretchinf function. Solving the roots of this equation gives the time stretching factor for non-linear file output cadence. + implicit none + ! Arguments + real(DP), intent(in) :: fstep_out + ! Result + real(DP) :: ans + + if (abs(fstep_out-1.0_DP) < epsilon(1.0_DP)) then + ans = self%nstep_out - tratio + else + ans = (1.0_DP - fstep_out**(self%nstep_out))/ (1.0_DP - fstep_out) - tratio + end if + + return + end function time_stretcher + end subroutine swiftest_io_param_reader From 620521a28545a81ecff097e35e7934920aa58c22 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 25 Feb 2023 20:11:29 -0500 Subject: [PATCH 3/3] Implemented time stretched output --- src/base/base_module.f90 | 1 + src/swiftest/swiftest_driver.f90 | 18 +++++++++++++++++- src/swiftest/swiftest_io.f90 | 10 ++++++++-- 3 files changed, 26 insertions(+), 3 deletions(-) diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 6060d53a1..480587b04 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -37,6 +37,7 @@ module base integer(I4B) :: istep_out = -1 !! Number of time steps between saved outputs integer(I4B) :: nstep_out = -1 !! Total number of saved outputs real(DP) :: fstep_out = 1.0_DP !! The output step time stretching factor + logical :: ltstretch = .false. !! Whether to employ time stretching or not character(STRMAX) :: outfile = BIN_OUTFILE !! Name of output binary file character(STRMAX) :: out_type = "NETCDF_DOUBLE" !! Binary format of output file character(STRMAX) :: out_form = "XVEL" !! Data to write to output file diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 8cc7f9a2a..2da264bea 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -27,6 +27,8 @@ program swiftest_driver integer(I8B) :: istart !! Starting index for loop counter integer(I4B) :: iout !! Output cadence counter integer(I4B) :: idump !! Dump cadence counter + integer(I4B) :: nout !! Current output step + integer(I4B) :: istep !! Current value of istep (used for time stretching) type(walltimer) :: integration_timer !! Object used for computing elapsed wall time call swiftest_io_get_args(integrator, param_file_name, display_style) @@ -44,6 +46,8 @@ program swiftest_driver iloop => param%iloop, & nloops => param%nloops, & istep_out => param%istep_out, & + fstep_out => param%fstep_out, & + ltstretch => param%ltstretch, & dump_cadence => param%dump_cadence, & display_unit => param%display_unit) @@ -53,6 +57,13 @@ program swiftest_driver iloop = istart - 1 iout = 0 idump = 0 + if (ltstretch) then + nout = floor(log(1.0_DP + iloop * (fstep_out - 1.0_DP) / istep_out) / log(fstep_out)) + istep = floor(istep_out * fstep_out**nout, kind=I4B) + else + nout = 1 + istep = istep_out + end if ! Set up nbody_system storage for intermittent file dumps if (dump_cadence == 0) dump_cadence = int(ceiling(nloops / (1.0_DP * istep_out), kind=I8B), kind=I4B) @@ -98,9 +109,14 @@ program swiftest_driver !> If the loop counter is at the output cadence value, append the data file with a single frame if (istep_out > 0) then iout = iout + 1 - if (iout == istep_out) then + if ((iout == istep) .or. (iloop == nloops)) then iout = 0 idump = idump + 1 + if (ltstretch) then + nout = nout + 1 + istep = floor(istep_out * fstep_out**nout, kind=I4B) + end if + call system_history%take_snapshot(param,nbody_system) if (idump == dump_cadence) then diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index b1af394ec..3e5efa606 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -2132,11 +2132,17 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i if (param%nstep_out <= 0) then param%nstep_out = int((param%tstop - param%t0) / (param%istep_out * param%dt)) param%fstep_out = 1.0_DP ! Linear output time + param%ltstretch = .false. else param%fstep_out = 1._DP tratio = (param%TSTOP - param%T0) / (param%istep_out * param%dt) - y = time_stretcher(param%fstep_out) - call solve_roots(time_stretcher,param%fstep_out) + if (int(tratio) == param%nstep_out) then + param%ltstretch = .false. + else + param%ltstretch = .true. + y = time_stretcher(param%fstep_out) + call solve_roots(time_stretcher,param%fstep_out) + end if end if if (param%dump_cadence < 0) then write(iomsg,*) 'Invalid DUMP_CADENCE. Must be a positive integer or 0.'