From 65eaabf8d698fcb3622f4e23c35db0f3fbf241cb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 25 Feb 2023 17:39:26 -0500 Subject: [PATCH] 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