Implemented time stretched output

MintonGroup · Feb 26, 2023 · 620521a · 620521a
1 parent 65eaabf
commit 620521a
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 3 deletions.
diff --git 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
@@ -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
@@ -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.'