From 9ce278fae18dc72bdd16fd623f30bb47935ea3c0 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Sun, 1 Mar 2020 17:18:56 +0100
Subject: [PATCH] [NL] INBTS; Add "fixed output times" option.

---
 .../t_fixed_output_times.md                   |  1 +
 ...CreateIterationNumberBasedTimeStepping.cpp | 16 +++++++-
 .../IterationNumberBasedTimeStepping.cpp      | 40 ++++++++++++++++++-
 .../IterationNumberBasedTimeStepping.h        | 11 ++++-
 .../TestTimeSteppingIterationNumber.cpp       | 33 ++++++++++++---
 5 files changed, 90 insertions(+), 11 deletions(-)
 create mode 100644 Documentation/ProjectFile/prj/time_loop/processes/process/time_stepping/IterationNumberBasedTimeStepping/t_fixed_output_times.md

diff --git a/Documentation/ProjectFile/prj/time_loop/processes/process/time_stepping/IterationNumberBasedTimeStepping/t_fixed_output_times.md b/Documentation/ProjectFile/prj/time_loop/processes/process/time_stepping/IterationNumberBasedTimeStepping/t_fixed_output_times.md
new file mode 100644
index 00000000000..c112eb6756d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/processes/process/time_stepping/IterationNumberBasedTimeStepping/t_fixed_output_times.md
@@ -0,0 +1 @@
+\copydoc NumLib::IterationNumberBasedTimeStepping::_fixed_output_times
diff --git a/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp
index 506b4c4199d..c57c891cecc 100644
--- a/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp
@@ -8,11 +8,12 @@
  */
 
 #include "CreateIterationNumberBasedTimeStepping.h"
+
 #include <string>
 
+#include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
-
 #include "IterationNumberBasedTimeStepping.h"
 #include "TimeStepAlgorithm.h"
 
@@ -43,8 +44,19 @@ std::unique_ptr<TimeStepAlgorithm> createIterationNumberBasedTimeStepping(
         //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__IterationNumberBasedTimeStepping__multiplier}
         config.getConfigParameter<std::vector<double>>("multiplier");
 
+    auto fixed_output_times =
+        //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__IterationNumberBasedTimeStepping__fixed_output_times}
+        config.getConfigParameter<std::vector<double>>("fixed_output_times",
+                                                       std::vector<double>{});
+    if (!fixed_output_times.empty())
+    {
+        // Remove possible duplicated elements and sort in descending order.
+        BaseLib::makeVectorUnique(fixed_output_times, std::greater<>());
+    }
+
     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),
+        std::move(fixed_output_times));
 }
 }  // namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
index 691b71b1647..497757b294c 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
@@ -18,16 +18,20 @@
 #include <limits>
 #include <utility>
 
+#include "BaseLib/Algorithm.h"
+
 namespace NumLib
 {
 IterationNumberBasedTimeStepping::IterationNumberBasedTimeStepping(
     double const t_initial, double const t_end, double const min_dt,
     double const max_dt, double const initial_dt,
     std::vector<int>&& iter_times_vector,
-    std::vector<double>&& multiplier_vector)
+    std::vector<double>&& multiplier_vector,
+    std::vector<double>&& fixed_output_times)
     : TimeStepAlgorithm(t_initial, t_end),
       _iter_times_vector(std::move(iter_times_vector)),
       _multiplier_vector(std::move(multiplier_vector)),
+      _fixed_output_times(std::move(fixed_output_times)),
       _min_dt(min_dt),
       _max_dt(max_dt),
       _initial_dt(initial_dt),
@@ -48,6 +52,9 @@ IterationNumberBasedTimeStepping::IterationNumberBasedTimeStepping(
     {
         OGS_FATAL("Vector of iteration numbers must be sorted.");
     }
+    //
+    // Remove possible duplicated elements and sort in descending order.
+    BaseLib::makeVectorUnique(_fixed_output_times, std::greater<double>());
 }
 
 bool IterationNumberBasedTimeStepping::next(double const /*solution_error*/,
@@ -72,7 +79,7 @@ bool IterationNumberBasedTimeStepping::next(double const /*solution_error*/,
 
     // prepare the next time step info
     _ts_current = _ts_prev;
-    _ts_current += getNextTimeStepSize();
+    _ts_current += checkSpecificTimeReached(getNextTimeStepSize());
 
     return true;
 }
@@ -126,6 +133,35 @@ double IterationNumberBasedTimeStepping::getNextTimeStepSize() const
     return dt;
 }
 
+double IterationNumberBasedTimeStepping::checkSpecificTimeReached(const double h_new)
+{
+    if (_fixed_output_times.empty())
+    {
+        return h_new;
+    }
+
+    const double specific_time = _fixed_output_times.back();
+    if ((specific_time > _ts_current.current()) &&
+        (_ts_current.current() + h_new - specific_time > 0.0))
+    {
+        _fixed_output_times.pop_back();
+        return specific_time - _ts_current.current();
+    }
+
+    return h_new;
+}
+
+void IterationNumberBasedTimeStepping::addFixedOutputTimes(
+    std::vector<double> const& extra_fixed_output_times)
+{
+    _fixed_output_times.insert(_fixed_output_times.end(),
+                               extra_fixed_output_times.begin(),
+                               extra_fixed_output_times.end());
+
+    // Remove possible duplicated elements and sort in descending order.
+    BaseLib::makeVectorUnique(_fixed_output_times, std::greater<double>());
+}
+
 bool IterationNumberBasedTimeStepping::accepted() const
 {
     return _accepted;
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h
index 78ef8a2c097..2bca95f977e 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h
@@ -81,6 +81,8 @@ public:
      * (\f$a_1\f$, \f$a_2\f$, ..., \f$a_n\f$) corresponding to the intervals
      * given by iter_times_vector.
      * A time step size is calculated by \f$\Delta t_{n+1} = a * \Delta t_{n}\f$
+     * @param fixed_output_times
+     * \copydoc NumLib::IterationNumberBasedTimeStepping::_fixed_output_times
      */
     IterationNumberBasedTimeStepping(double const t_initial,
                                      double const t_end,
@@ -88,7 +90,8 @@ public:
                                      double const max_dt,
                                      double const initial_dt,
                                      std::vector<int>&& iter_times_vector,
-                                     std::vector<double>&& multiplier_vector);
+                                     std::vector<double>&& multiplier_vector,
+                                     std::vector<double>&& fixed_output_times);
 
     ~IterationNumberBasedTimeStepping() override = default;
 
@@ -104,6 +107,8 @@ public:
     /// Return the number of repeated steps.
     int getNumberOfRepeatedSteps() const { return _n_rejected_steps; }
 
+    void addFixedOutputTimes(
+        std::vector<double> const& extra_fixed_output_times) override;
 private:
     /// Calculate the next time step size.
     double getNextTimeStepSize() const;
@@ -111,11 +116,15 @@ private:
     /// Find a multiplier for the given number of iterations.
     double findMultiplier(int number_iterations) const;
 
+    double checkSpecificTimeReached(const double h_new);
+
     /// This vector stores the number of iterations to which the respective
     /// multiplier coefficient will be applied.
     const std::vector<int> _iter_times_vector;
     /// This vector stores the multiplier coefficients.
     const std::vector<double> _multiplier_vector;
+    /// Timesteps to be taken independent of the time stepping scheme.
+    std::vector<double> _fixed_output_times;
     /// The minimum allowed time step size.
     const double _min_dt;
     /// The maximum allowed time step size.
diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
index 24a83dd91f8..d9c9ab49fab 100644
--- a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
+++ b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
@@ -26,9 +26,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
 {
     std::vector<int> iter_times_vector = {0, 3, 5, 7};
     std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
-    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1,
-                                                 std::move(iter_times_vector),
-                                                 std::move(multiplier_vector));
+    NumLib::IterationNumberBasedTimeStepping alg(
+        1, 31, 1, 10, 1, std::move(iter_times_vector),
+        std::move(multiplier_vector), {});
 
     const double solution_error = 0.;
 
@@ -87,9 +87,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased2)
 {
     std::vector<int> iter_times_vector = {0, 3, 5, 7};
     std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
-    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1,
-                                                 std::move(iter_times_vector),
-                                                 std::move(multiplier_vector));
+    NumLib::IterationNumberBasedTimeStepping alg(
+        1, 31, 1, 10, 1, std::move(iter_times_vector),
+        std::move(multiplier_vector), {});
 
     std::vector<int> nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1};
     const std::vector<double> expected_vec_t = {1,  2,  4,  8,  16,
@@ -101,3 +101,24 @@ TEST(NumLib, TimeSteppingIterationNumberBased2)
     ASSERT_EQ(0u, alg.getNumberOfRepeatedSteps());
     ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
 }
+
+TEST(NumLib, TimeSteppingIterationNumberBased2FixedOutputTimes)
+{
+    std::vector<int> iter_times_vector = {0, 3, 5, 7};
+    std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
+    std::vector<double> fixed_output_times = {5, 20};
+    NumLib::IterationNumberBasedTimeStepping alg(
+        1, 31, 1, 10, 1, std::move(iter_times_vector),
+        std::move(multiplier_vector), std::move(fixed_output_times));
+
+    std::vector<int> nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1, 1, 1, 1, 1};
+    const std::vector<double> expected_vec_t = {1,  2,  4,  5,  7,  9,  10,
+                                                11, 12, 14, 18, 20, 24, 31};
+
+    std::vector<double> vec_t = timeStepping(alg, nr_iterations);
+
+    EXPECT_EQ(expected_vec_t.size(), vec_t.size());
+    ASSERT_EQ(0u, alg.getNumberOfRepeatedSteps());
+    ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                      std::numeric_limits<double>::epsilon());
+}
-- 
GitLab