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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Feb 26, 2023
2 parents c34cecb + 620521a commit e26b277
Show file tree
Hide file tree
Showing 6 changed files with 297 additions and 55 deletions.
34 changes: 29 additions & 5 deletions python/swiftest/swiftest/simulation_class.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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",
}

Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand Down
3 changes: 3 additions & 0 deletions src/base/base_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,9 @@ 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
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
Expand Down
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.1_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 * (1.0_DP - 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 = 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
Loading

0 comments on commit e26b277

Please sign in to comment.