[NL] Include fixed output times in FixedTimeStepping algorithm

NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp

@@ -11,8 +11,13 @@
 
 #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"
@@ -20,6 +25,102 @@
 
 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,
@@ -33,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");
@@ -64,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)
             {
@@ -100,10 +201,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)
         {
@@ -128,6 +229,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

@@ -23,6 +23,10 @@ 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(