Merge pull request #2831 from endJunction/FixedOutputTimesForIterationNumberBasedTS

Fixed output times for iteration number based time stepping algorithm

Documentation/ProjectFile/prj/time_loop/processes/process/time_stepping/IterationNumberBasedTimeStepping/t_fixed_output_times.md

+\copydoc NumLib::IterationNumberBasedTimeStepping::_fixed_output_times

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

NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp

@@ -20,6 +20,26 @@
 
 namespace NumLib
 {
+EvolutionaryPIDcontroller::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,
+    std::vector<double>&& fixed_output_times, const double tol)
+    : TimeStepAlgorithm(t0, t_end),
+      _h0(h0),
+      _h_min(h_min),
+      _h_max(h_max),
+      _rel_h_min(rel_h_min),
+      _rel_h_max(rel_h_max),
+      _fixed_output_times(std::move(fixed_output_times)),
+      _tol(tol),
+      _e_n_minus1(0.),
+      _e_n_minus2(0.),
+      _is_accepted(true)
+{
+    // Remove possible duplicated elements. Result will be sorted.
+    BaseLib::makeVectorUnique(_fixed_output_times);
+}
+
 bool EvolutionaryPIDcontroller::next(double const solution_error,
                                      int const /*number_iterations*/)
 {
@@ -36,7 +56,8 @@ bool EvolutionaryPIDcontroller::next(double const solution_error,
                                               : 0.5 * _ts_current.dt();
 
         h_new = limitStepSize(h_new, is_previous_step_accepted);
-        h_new = checkSpecificTimeReached(h_new);
+        h_new = possiblyClampDtToNextFixedTime(_ts_current.current(), h_new,
+                                               _fixed_output_times);
 
         _ts_current = _ts_prev;
         _ts_current += h_new;
@@ -96,7 +117,8 @@ bool EvolutionaryPIDcontroller::next(double const solution_error,
         }
 
         h_new = limitStepSize(h_new, is_previous_step_accepted);
-        h_new = checkSpecificTimeReached(h_new);
+        h_new = possiblyClampDtToNextFixedTime(_ts_current.current(), h_new,
+                                               _fixed_output_times);
         _dt_vector.push_back(h_new);
 
         _ts_prev = _ts_current;
@@ -145,24 +167,6 @@ double EvolutionaryPIDcontroller::limitStepSize(
     return limited_h;
 }
 
-double EvolutionaryPIDcontroller::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 EvolutionaryPIDcontroller::addFixedOutputTimes(
     std::vector<double> const& extra_fixed_output_times)
 {
@@ -170,8 +174,8 @@ void EvolutionaryPIDcontroller::addFixedOutputTimes(
                                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>());
+    // Remove possible duplicated elements. Result will be sorted.
+    BaseLib::makeVectorUnique(_fixed_output_times);
 }
 
 bool EvolutionaryPIDcontroller::canReduceTimestepSize() const

NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h

@@ -55,21 +55,8 @@ public:
                               const double h0, const double h_min,
                               const double h_max, const double rel_h_min,
                               const double rel_h_max,
-                              const std::vector<double>&& fixed_output_times,
-                              const double tol)
-        : TimeStepAlgorithm(t0, t_end),
-          _h0(h0),
-          _h_min(h_min),
-          _h_max(h_max),
-          _rel_h_min(rel_h_min),
-          _rel_h_max(rel_h_max),
-          _fixed_output_times(std::move(fixed_output_times)),
-          _tol(tol),
-          _e_n_minus1(0.),
-          _e_n_minus2(0.),
-          _is_accepted(true)
-    {
-    }
+                              std::vector<double>&& fixed_output_times,
+                              const double tol);
 
     bool next(double solution_error, int number_iterations) override;
 
@@ -124,8 +111,6 @@ private:
      */
     double limitStepSize(const double h_new,
                          const bool previous_step_accepted) const;
-
-    double checkSpecificTimeReached(const double h_new);
 };
 
 /// Create an EvolutionaryPIDcontroller time stepper from the given

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. Result will be sorted.
+    BaseLib::makeVectorUnique(_fixed_output_times);
 }
 
 bool IterationNumberBasedTimeStepping::next(double const /*solution_error*/,
@@ -72,7 +79,8 @@ bool IterationNumberBasedTimeStepping::next(double const /*solution_error*/,
 
     // prepare the next time step info
     _ts_current = _ts_prev;
-    _ts_current += getNextTimeStepSize();
+    _ts_current += possiblyClampDtToNextFixedTime(
+        _ts_current.current(), getNextTimeStepSize(), _fixed_output_times);
 
     return true;
 }
@@ -126,6 +134,17 @@ double IterationNumberBasedTimeStepping::getNextTimeStepSize() const
     return dt;
 }
 
+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. Result will be sorted.
+    BaseLib::makeVectorUnique(_fixed_output_times);
+}
+
 bool IterationNumberBasedTimeStepping::accepted() const
 {
     return _accepted;

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;
@@ -116,6 +121,8 @@ private:
     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.

NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "TimeStepAlgorithm.h"
+
+#include <algorithm>
+
+namespace NumLib
+{
+double possiblyClampDtToNextFixedTime(
+    double const t, double const dt,
+    std::vector<double> const& fixed_output_times)
+{
+    auto const specific_time = std::upper_bound(
+        std::cbegin(fixed_output_times), std::cend(fixed_output_times), t);
+
+    if (specific_time == std::cend(fixed_output_times))
+    {
+        return dt;
+    }
+
+    if ((*specific_time > t) && (t + dt - *specific_time > 0.0))
+    {
+        return *specific_time - t;
+    }
+
+    return dt;
+}
+}  // namespace NumLib

NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h

@@ -141,4 +141,14 @@ protected:
     std::vector<double> _dt_vector;
 };
 
+/// If any of the fixed times will be reached with given time increment, it will
+/// be reduced, otherwise the input will be returned.
+/// \pre The input vector of fixed times must be sorted.
+/// \param t Current time.
+/// \param dt Suggested time increment.
+/// \param fixed_output_times Sorted list of times which are to be reached.
+double possiblyClampDtToNextFixedTime(
+    double const t, double const dt,
+    std::vector<double> const& fixed_output_times);
+
 }  // namespace NumLib

Tests/Data/Parabolic/ComponentTransport/MassConservation/mass_conservation.prj

@@ -124,7 +124,6 @@
                     <each_steps>300</each_steps>
                 </pair>
             </timesteps>
-            <fixed_output_times> 434629.8 866629.8 1298629.8 1730629.8 2162629.8 2594629.8 </fixed_output_times>
             <output_iteration_results>false</output_iteration_results>
             <variables>
                 <variable>concentration</variable>

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());
+}