Merge branch 'PassFixedOutputTimesToTimeStepper' into 'master'

Pass fixed output times to time stepper

See merge request ogs/ogs!4652

NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp

@@ -20,7 +20,8 @@ namespace NumLib
 {
 class TimeStepAlgorithm;
 std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(
-    BaseLib::ConfigTree const& config)
+    BaseLib::ConfigTree const& config,
+    std::vector<double> const& fixed_times_for_output)
 {
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__type}
     config.checkConfigParameter("type", "EvolutionaryPIDcontroller");
@@ -44,7 +45,14 @@ std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__EvolutionaryPIDcontroller__tol}
     auto const tol = config.getConfigParameter<double>("tol");
 
-    return std::make_unique<EvolutionaryPIDcontroller>(
-        t0, t_end, h0, h_min, h_max, rel_h_min, rel_h_max, tol);
+    return std::make_unique<EvolutionaryPIDcontroller>(t0,
+                                                       t_end,
+                                                       h0,
+                                                       h_min,
+                                                       h_max,
+                                                       rel_h_min,
+                                                       rel_h_max,
+                                                       tol,
+                                                       fixed_times_for_output);
 }
 }  // end of namespace NumLib

NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h

@@ -12,6 +12,7 @@
 #pragma once
 
 #include <memory>
+#include <vector>
 
 namespace BaseLib
 {
@@ -25,5 +26,6 @@ class TimeStepAlgorithm;
 /// Create an EvolutionaryPIDcontroller time stepper from the given
 /// configuration
 std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(
-    BaseLib::ConfigTree const& config);
+    BaseLib::ConfigTree const& config,
+    std::vector<double> const& fixed_times_for_output);
 }  // end of namespace NumLib

NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp

@@ -11,18 +11,120 @@
 
 #include "CreateFixedTimeStepping.h"
 
+#include <fmt/ranges.h>
+
+#include <algorithm>
+#include <numeric>
 #include <string>
 
+#include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
 #include "FixedTimeStepping.h"
 #include "TimeStepAlgorithm.h"
 
 namespace NumLib
 {
+std::size_t findDeltatInterval(double const t_initial,
+                               std::vector<double> const& delta_ts,
+                               double const fixed_output_time)
+{
+    if (fixed_output_time < t_initial)
+    {
+        return std::numeric_limits<std::size_t>::max();
+    }
+
+    auto timestepper_time = t_initial;
+    for (std::size_t k = 0; k < delta_ts.size(); ++k)
+    {
+        if (timestepper_time <= fixed_output_time &&
+            fixed_output_time < timestepper_time + delta_ts[k])
+        {
+            return k;
+        }
+        timestepper_time += delta_ts[k];
+    }
+    if (fixed_output_time == timestepper_time + delta_ts.back())
+    {
+        return std::numeric_limits<std::size_t>::max();
+    }
+    return std::numeric_limits<std::size_t>::max();
+}
+
+void incorporateFixedTimesForOutput(
+    double const t_initial, double const t_end, std::vector<double>& delta_ts,
+    std::vector<double> const& fixed_times_for_output)
+{
+    if (fixed_times_for_output.empty())
+    {
+        return;
+    }
+
+    if (auto lower_bound =
+            std::lower_bound(begin(fixed_times_for_output),
+                             end(fixed_times_for_output), t_initial);
+        lower_bound != begin(fixed_times_for_output))
+    {
+        WARN(
+            "Request for output at times {}, but the simulation's start time "
+            "is {}. Output will be skipped.",
+            fmt::join(begin(fixed_times_for_output), lower_bound, ", "),
+            t_initial);
+    }
+
+    if (auto upper_bound = std::upper_bound(begin(fixed_times_for_output),
+                                            end(fixed_times_for_output), t_end);
+        upper_bound != end(fixed_times_for_output))
+    {
+        WARN(
+            "Request for output at times {}, but simulation's end time is {}. "
+            "Output will be skipped.",
+            fmt::join(upper_bound, end(fixed_times_for_output), ", "),
+            t_end);
+    }
+
+    if (delta_ts.empty())
+    {
+        WARN("No timesteps specified.");
+        return;
+    }
+
+    // incorporate fixed output times into dts vector
+    for (auto const fixed_time_for_output : fixed_times_for_output)
+    {
+        auto const interval_number =
+            findDeltatInterval(t_initial, delta_ts, fixed_time_for_output);
+        if (interval_number == std::numeric_limits<std::size_t>::max())
+        {
+            WARN("Did not find interval for fixed output time {}",
+                 fixed_time_for_output);
+            continue;
+        }
+        auto const lower_bound = std::accumulate(
+            begin(delta_ts), begin(delta_ts) + interval_number, t_initial);
+        auto const upper_bound = std::accumulate(
+            begin(delta_ts), begin(delta_ts) + interval_number + 1, t_initial);
+        if (fixed_time_for_output - lower_bound <=
+            std::numeric_limits<double>::epsilon())
+        {
+            continue;
+        }
+        if (upper_bound - fixed_time_for_output <=
+            std::numeric_limits<double>::epsilon())
+        {
+            continue;
+        }
+        delta_ts[interval_number] = fixed_time_for_output - lower_bound;
+
+        delta_ts.insert(delta_ts.begin() + interval_number + 1,
+                        upper_bound - fixed_time_for_output);
+    }
+}
+
 class TimeStepAlgorithm;
 std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
-    BaseLib::ConfigTree const& config)
+    BaseLib::ConfigTree const& config,
+    std::vector<double> const& fixed_times_for_output)
 {
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__type}
     config.checkConfigParameter("type", "FixedTimeStepping");
@@ -32,15 +134,15 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__FixedTimeStepping__t_end}
     auto const t_end = config.getConfigParameter<double>("t_end");
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__FixedTimeStepping__timesteps}
-    auto const delta_ts = config.getConfigSubtree("timesteps");
+    auto const delta_ts_config = config.getConfigSubtree("timesteps");
 
-    std::vector<double> timesteps;
+    std::vector<double> delta_ts;
     double t_curr = t_initial;
     double delta_t = 0.0;
 
     // TODO: consider adding call "listNonEmpty" to config tree
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__FixedTimeStepping__timesteps__pair}
-    auto const range = delta_ts.getConfigSubtreeList("pair");
+    auto const range = delta_ts_config.getConfigSubtreeList("pair");
     if (range.begin() == range.end())
     {
         OGS_FATAL("no timesteps have been given");
@@ -63,10 +165,10 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
 
         if (t_curr <= t_end)
         {
-            auto const new_size = timesteps.size() + repeat;
+            auto const new_size = delta_ts.size() + repeat;
             try
             {
-                timesteps.resize(new_size, delta_t);
+                delta_ts.resize(new_size, delta_t);
             }
             catch (std::length_error const& e)
             {
@@ -90,7 +192,11 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
                     e.what());
             }
 
-            t_curr += repeat * delta_t;
+            // Multiplying dt * repeat is not the same as in the current
+            // implementation of the time loop, where the dt's are added.
+            // Therefore the sum of all dt is taken here.
+            t_curr =
+                std::accumulate(end(delta_ts) - repeat, end(delta_ts), t_curr);
         }
     }
 
@@ -99,10 +205,10 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
     {
         auto const repeat =
             static_cast<std::size_t>(std::ceil((t_end - t_curr) / delta_t));
-        auto const new_size = timesteps.size() + repeat;
+        auto const new_size = delta_ts.size() + repeat;
         try
         {
-            timesteps.resize(new_size, delta_t);
+            delta_ts.resize(new_size, delta_t);
         }
         catch (std::length_error const& e)
         {
@@ -127,6 +233,8 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
         }
     }
 
-    return std::make_unique<FixedTimeStepping>(t_initial, t_end, timesteps);
+    incorporateFixedTimesForOutput(t_initial, t_end, delta_ts,
+                                   fixed_times_for_output);
+    return std::make_unique<FixedTimeStepping>(t_initial, t_end, delta_ts);
 }
 }  // end of namespace NumLib

NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h

@@ -12,6 +12,7 @@
 #pragma once
 
 #include <memory>
+#include <vector>
 
 namespace BaseLib
 {
@@ -22,8 +23,13 @@ namespace NumLib
 {
 class TimeStepAlgorithm;
 
+std::size_t findDeltatInterval(double const t_initial,
+                               std::vector<double> const& delta_ts,
+                               double const fixed_output_time);
+
 /// Create a FixedTimeStepping time stepper from the given
 /// configuration
 std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
-    BaseLib::ConfigTree const& config);
+    BaseLib::ConfigTree const& config,
+    std::vector<double> const& fixed_times_for_output);
 }  // end of namespace NumLib

NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp

@@ -20,7 +20,8 @@ namespace NumLib
 {
 class TimeStepAlgorithm;
 std::unique_ptr<TimeStepAlgorithm> createIterationNumberBasedTimeStepping(
-    BaseLib::ConfigTree const& config)
+    BaseLib::ConfigTree const& config,
+    std::vector<double> const& fixed_times_for_output)
 {
     //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__type}
     config.checkConfigParameter("type", "IterationNumberBasedTimeStepping");
@@ -45,6 +46,7 @@ std::unique_ptr<TimeStepAlgorithm> createIterationNumberBasedTimeStepping(
 
     return std::make_unique<IterationNumberBasedTimeStepping>(
         t_initial, t_end, minimum_dt, maximum_dt, initial_dt,
-        std::move(number_iterations), std::move(multiplier));
+        std::move(number_iterations), std::move(multiplier),
+        fixed_times_for_output);
 }
 }  // namespace NumLib

NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.h

@@ -10,6 +10,7 @@
 #pragma once
 
 #include <memory>
+#include <vector>
 
 namespace BaseLib
 {
@@ -26,5 +27,6 @@ namespace NumLib
 /// Create a IterationNumberBasedTimeStepping time stepper from the given
 /// configuration.
 std::unique_ptr<TimeStepAlgorithm> createIterationNumberBasedTimeStepping(
-    BaseLib::ConfigTree const& config);
+    BaseLib::ConfigTree const& config,
+    std::vector<double> const& fixed_times_for_output);
 }  // namespace NumLib

NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp

@@ -136,6 +136,34 @@ double EvolutionaryPIDcontroller::limitStepSize(
             limited_h = std::max(_h_min, 0.5 * timestep_current.dt());
         }
     }
+
+    if (_fixed_times_for_output.empty())
+    {
+        return limited_h;
+    }
+
+    // find first fixed timestep for output larger than the current time, i.e.,
+    // current time < fixed output time
+    auto fixed_output_time_it = std::find_if(
+        std::begin(_fixed_times_for_output), std::end(_fixed_times_for_output),
+        [&timestep_current](auto const fixed_output_time)
+        { return timestep_current.current() < fixed_output_time; });
+
+    if (fixed_output_time_it != _fixed_times_for_output.end())
+    {
+        // check if the fixed output time is in the interval
+        // (current time, current time + limited_h)
+        if (*fixed_output_time_it < timestep_current.current() + limited_h)
+        {
+            // check if the potential adjusted time step is larger than zero
+            if (std::abs(*fixed_output_time_it - timestep_current.current()) >
+                std::numeric_limits<double>::epsilon() *
+                    timestep_current.current())
+            {
+                return *fixed_output_time_it - timestep_current.current();
+            }
+        }
+    }
     return limited_h;
 }
 

NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h

@@ -49,7 +49,8 @@ class EvolutionaryPIDcontroller final : public TimeStepAlgorithm
     EvolutionaryPIDcontroller(const double t0, const double t_end,
                               const double h0, const double h_min,
                               const double h_max, const double rel_h_min,
-                              const double rel_h_max, const double tol)
+                              const double rel_h_max, const double tol,
+                              std::vector<double> const& fixed_times_for_output)
         : TimeStepAlgorithm(t0, t_end),
           _h0(h0),
           _h_min(h_min),
@@ -58,7 +59,8 @@ class EvolutionaryPIDcontroller final : public TimeStepAlgorithm
           _rel_h_max(rel_h_max),
           _tol(tol),
           _e_n_minus1(0.),
-          _e_n_minus2(0.)
+          _e_n_minus2(0.),
+          _fixed_times_for_output(fixed_times_for_output)
     {
     }
 
@@ -92,6 +94,8 @@ class EvolutionaryPIDcontroller final : public TimeStepAlgorithm
     double _e_n_minus1;  ///< \f$e_{n-1}\f$.
     double _e_n_minus2;  ///< \f$e_{n-2}\f$.
 
+    std::vector<double> const _fixed_times_for_output;
+
     /**
      * Force the computed time step size in the given range
      * (see the formulas in the documentation of the class)