Merge branch 'IterationNumberBasedTimeStepping' into 'master'

Fix bugs and add tests for IterationNumberBasedTimeStepping

See merge request ogs/ogs!5383

Author: TobiasMeisel

NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp

@@ -95,6 +95,8 @@ std::tuple<bool, double> IterationNumberBasedTimeStepping::next(
         ts_previous =  // essentially equal to ts_previous.dt = _ts_current.dt.
             TimeStep{ts_previous.previous(), ts_previous.previous() + dt,
                      ts_previous.timeStepNumber()};
+        ts_current = TimeStep{ts_current.previous(), ts_current.previous() + dt,
+                              ts_current.timeStepNumber()};
 
         _previous_time_step_accepted = false;
         return std::make_tuple(_previous_time_step_accepted, dt);
@@ -165,6 +167,9 @@ double IterationNumberBasedTimeStepping::getNextTimeStepSize(
         return std::clamp(dt, _min_dt, _max_dt);
     }
 
+    // restrict dt to _max_dt before taking fixed times for output into account
+    dt = std::min(dt, _max_dt);
+
     // 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(

NumLib/TimeStepping/TimeStep.h

@@ -12,6 +12,8 @@
 
 #pragma once
 
+#include <spdlog/fmt/bundled/ostream.h>
+
 #include <cstddef>
 #include <limits>
 
@@ -84,6 +86,14 @@ class TimeStep final
     void setAccepted(bool const accepted) { _is_accepted = accepted; }
     bool isAccepted() const { return _is_accepted; }
 
+    friend inline std::ostream& operator<<(std::ostream& os, TimeStep const& ts)
+    {
+        return os << "previous: " << ts.previous()
+                  << " | current: " << ts.current() << " | dt: " << ts.dt()
+                  << " | timestep number: " << ts.timeStepNumber()
+                  << " | is_accepted: " << ts.isAccepted();
+    }
+
 private:
     /// previous time step
     Time _previous;
@@ -105,3 +115,8 @@ inline void updateTimeSteps(double const dt, TimeStep& previous_timestep,
 }
 
 }  // namespace NumLib
+
+template <>
+struct fmt::formatter<NumLib::TimeStep> : fmt::ostream_formatter
+{
+};

ProcessLib/SmallDeformation/Tests.cmake

@@ -363,9 +363,7 @@ AddTest(
     REQUIREMENTS NOT OGS_USE_MPI
     RUNTIME 1
     PROPERTIES
-    # Time step size of 0.125 means the time step has been halved twice from its initial size.
-    # That proves that OGS gracefully handles such conditions
-    PASS_REGULAR_EXPRESSION "The new step size of 0[.]125 is the same as that of the previous rejected time step[.]"
+    PASS_REGULAR_EXPRESSION "Time stepper cannot reduce the time step size further[.]"
 )
 
 if(NOT OGS_USE_PETSC)

Tests/Data/TH2M/H2M/EmbeddedFracturePermeability/IfG_t_0.vtu

@@ -83,7 +83,7 @@ class NumLibTimeSteppingFixed_TimeSteps : public ::testing::Test
                             std::vector<double> const expected_times) const
     {
         std::vector<double> times =
-            timeStepping(algorithm, dummy_number_iterations, {});
+            timeStepping(algorithm, dummy_number_iterations, {}, {});
 
         ASSERT_EQ(expected_times.size(), times.size());
         ASSERT_ARRAY_NEAR(expected_times, times, expected_times.size(),

Tests/Data/TH2M/H2M/EmbeddedFracturePermeability/IfG_t_100.vtu

@@ -12,7 +12,11 @@
 
 #include <gtest/gtest.h>
 
-#include <ranges>
+#include <algorithm>
+#include <range/v3/range/conversion.hpp>
+#include <range/v3/view/iota.hpp>
+#include <range/v3/view/transform.hpp>
+#include <range/v3/view/zip.hpp>
 #include <tuple>
 #include <utility>
 #include <vector>
@@ -71,7 +75,7 @@ TEST(NumLib, findMultiplierTimestepAcceptedPiecewiseConstant)
 
     bool const current_time_step_is_accepted = true;
     for (auto const [n, m] :
-         std::views::zip(nonlinear_iteration_numbers, multipliers))
+         ranges::views::zip(nonlinear_iteration_numbers, multipliers))
     {
         auto const multiplier = NumLib::findMultiplier(
             n, current_time_step_is_accepted, nonlinear_iteration_numbers,
@@ -91,7 +95,7 @@ TEST(
 
     bool const current_time_step_is_accepted = true;
     for (auto const [n, m] :
-         std::views::zip(nonlinear_iteration_numbers, multipliers))
+         ranges::views::zip(nonlinear_iteration_numbers, multipliers))
     {
         auto const multiplier = NumLib::findMultiplier(
             n, current_time_step_is_accepted, nonlinear_iteration_numbers,
@@ -187,7 +191,7 @@ TEST(NumLib, findMultiplierTimestepRejected)
 
     bool const current_time_step_is_accepted = false;
     for (auto const [n, m] :
-         std::views::zip(nonlinear_iteration_numbers, multipliers))
+         ranges::views::zip(nonlinear_iteration_numbers, multipliers))
     {
         auto const multiplier = NumLib::findMultiplier(
             n, current_time_step_is_accepted, nonlinear_iteration_numbers,
@@ -328,11 +332,11 @@ TEST(NumLib, TimeSteppingIterationNumberBased2)
         1, 31, 1, 10, 1, multiplier_interpolation_type,
         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,
+    std::vector<int> const nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1};
+    std::vector<double> const expected_vec_t = {1,  2,  4,  8,  16,
                                                 24, 28, 29, 30, 31};
 
-    std::vector<double> vec_t = timeStepping(alg, nr_iterations, {});
+    std::vector<double> vec_t = timeStepping(alg, nr_iterations, {}, {});
 
     ASSERT_EQ(expected_vec_t.size(), vec_t.size());
     ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
@@ -350,14 +354,196 @@ TEST(NumLib, TimeSteppingIterationNumberBased2FixedOutputTimes)
         1, 31, 1, 10, 1, multiplier_interpolation_type,
         std::move(iter_times_vector), std::move(multiplier_vector), {});
 
-    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,
+    std::vector<int> const nr_iterations = {0, 2, 2, 2, 4, 6, 8,
+                                            4, 1, 1, 1, 1, 1};
+    std::vector<double> const 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, fixed_output_times);
+        timeStepping(alg, nr_iterations, fixed_output_times, {});
 
     EXPECT_EQ(expected_vec_t.size(), vec_t.size());
     ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
                       std::numeric_limits<double>::epsilon());
 }
+
+TEST(NumLib, TimeSteppingIterationNumberBased_simple)
+{
+    // *** initialization of IterationNumberBaseTimeStepping object
+    constexpr int number_of_multipliers = 20;
+    std::vector multiplier_vector(number_of_multipliers, 1.0);
+
+    auto iter_times_vector =
+        ranges::views::iota(1, static_cast<int>(multiplier_vector.size() + 1)) |
+        ranges::to<std::vector>;
+
+    std::vector<double> const fixed_output_times = {};
+
+    NumLib::MultiplyerInterpolationType const multiplier_interpolation_type =
+        NumLib::MultiplyerInterpolationType::PiecewiseConstant;
+
+    double const t_initial = 0.0;
+    double const t_end = 8000.0;
+    double const min_dt = 0.001;
+    double const max_dt = 2100;
+    double const initial_dt = 100;
+
+    NumLib::IterationNumberBasedTimeStepping alg(
+        t_initial, t_end, min_dt, max_dt, initial_dt,
+        multiplier_interpolation_type, std::move(iter_times_vector),
+        std::move(multiplier_vector), std::move(fixed_output_times));
+    // *** end initialization of IterationNumberBaseTimeStepping object
+
+    std::vector<int> const rejected_steps = {};
+    std::vector<int> const nr_iterations = {
+        0, 4, 3,  9, 3,  5, 4, 3, 3, 7, 10, 4, 3, 3, 3, 5, 3, 4, 7, 15,
+        3, 7, 15, 3, 11, 8, 4, 5, 5, 4, 3,  3, 3, 3, 3, 3, 3, 5, 3, 3,
+        3, 3, 3,  3, 7,  5, 4, 4, 3, 3, 3,  3, 3, 3, 3, 3, 3, 3, 3, 4,
+        4, 3, 3,  3, 3,  3, 3, 3, 3, 3, 3,  3, 3, 3, 3, 3, 4, 3, 3, 3};
+
+    auto const expected_vec_t =
+        ranges::views::iota(0, static_cast<int>(nr_iterations.size()) + 1) |
+        ranges::views::transform([](int i) { return (i * 100.0); }) |
+        ranges::to<std::vector>;
+
+    std::vector<double> const vec_t =
+        timeStepping(alg, nr_iterations, fixed_output_times, rejected_steps);
+
+    EXPECT_EQ(expected_vec_t.size(), vec_t.size());
+    ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                      std::numeric_limits<double>::epsilon());
+}
+
+TEST(NumLib, TimeSteppingIterationNumberBased_simple2)
+{
+    // *** initialization of IterationNumberBaseTimeStepping object
+    constexpr int number_of_multipliers = 20;
+    std::vector multiplier_vector(number_of_multipliers - 1, 1.0);
+    multiplier_vector.emplace_back(0.1);  // multiplier for rejected step
+
+    auto iter_times_vector =
+        ranges::views::iota(1, static_cast<int>(multiplier_vector.size() + 1)) |
+        ranges::to<std::vector>;
+
+    std::vector<double> fixed_output_times = {};
+
+    NumLib::MultiplyerInterpolationType const multiplier_interpolation_type =
+        NumLib::MultiplyerInterpolationType::PiecewiseConstant;
+    double const t_initial = 0.0;
+    double const t_end = 200.0;
+    double const min_dt = 10;
+    double const max_dt = 200;
+    double const initial_dt = 100;
+    NumLib::IterationNumberBasedTimeStepping alg(
+        t_initial, t_end, min_dt, max_dt, initial_dt,
+        multiplier_interpolation_type, std::move(iter_times_vector),
+        std::move(multiplier_vector), std::move(fixed_output_times));
+    // *** end initialization of IterationNumberBaseTimeStepping object
+
+    std::vector<int> const rejected_steps = {2};
+    std::vector<int> const nr_iterations = {0, 4, 3, 3,  3,  3,
+                                            3, 3, 9, 10, 19, 1};
+    // current time step size:
+    std::vector<double> const expected_vec_t = {
+        0, 100, 200, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200};
+
+    std::vector<double> const vec_t =
+        timeStepping(alg, nr_iterations, fixed_output_times, rejected_steps);
+
+    EXPECT_EQ(expected_vec_t.size(), vec_t.size());
+    ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                      std::numeric_limits<double>::epsilon());
+}
+
+// extracted from Sandwich-0c-Aqeel-R5DeCoT68
+TEST(NumLib, TimeSteppingIterationNumberBasedSandwich_0c)
+{
+    // *** initialization of IterationNumberBaseTimeStepping object
+    std::vector<double> multiplier_vector = {
+        5,      4.25,   3.5,    2.75,   2,      1.6333, 1.2667,
+        0.9,    0.8333, 0.7667, 0.7,    0.6333, 0.5667, 0.5,
+        0.4333, 0.3667, 0.3,    0.2333, 0.1667, 0.1};
+
+    auto iter_times_vector =
+        ranges::views::iota(1, static_cast<int>(multiplier_vector.size() + 1)) |
+        ranges::to<std::vector>;
+
+    std::vector<double> fixed_output_times = {
+        0,       432000,  440640,  950400,  959040,  3024000, 3032640,
+        5961600, 6480000, 7689600, 7776000, 7819200, 12528000};
+
+    NumLib::MultiplyerInterpolationType const multiplier_interpolation_type =
+        NumLib::MultiplyerInterpolationType::PiecewiseConstant;
+    // *** original data
+    double const t_initial = 0.0;
+    double const t_end = 12528000.0;
+    double const min_dt = 0.001;
+    double const max_dt = 432000;
+    double const initial_dt = 120;
+    NumLib::IterationNumberBasedTimeStepping alg(
+        t_initial, t_end, min_dt, max_dt, initial_dt,
+        multiplier_interpolation_type, std::move(iter_times_vector),
+        std::move(multiplier_vector), std::move(fixed_output_times));
+    // *** end initialization of IterationNumberBaseTimeStepping object
+
+    std::vector<int> const rejected_steps = {6, 8, 10, 12, 14, 14, 24};
+
+    std::vector<int> const nr_iterations = {
+        0,  4, 3, 9, 3, 5, 5,  // <- time step 6 rejected
+        3,                     // <- time step 6 repeated
+        3,  8,                 // <- time step 8 rejected
+        10,                    // <- time step 8 repeated
+        4,  4,                 // <- time step 10 rejected
+        3,                     // <- time step 10 repeated
+        3,  6,                 // <- time step 12 rejected
+        3,                     // <- time step 12 repeated
+        4,  8,                 // <- time step 14 rejected
+        17,                    // <- time step 14 rejected again
+        3,                     // <- time step 14 repeated (and accepted)
+        9,  4, 5, 5, 4, 3, 3, 3, 3, 11,  // <- time step 24 rejected
+        3                                // <- time step 24 repeated
+    };
+
+    // current timestep sizes:
+    std::vector<double> const expected_vec_t = {
+        0,
+        120,
+        450,
+        1605,
+        2567.4614999999999,
+        5936.0767500000002,
+        12673.30725,         // ts 6 rejected
+        6609.7998000000007,  // ts 6 repeated
+        8967.8304749999988,
+        17220.937837499994,  // ts 8 rejected
+        16395.627101249993,  // ts 8 repeated
+        22090.518774595865,
+        37751.470876297011,  // ts 10 rejected
+        23656.61398476598,   // ts 10 repeated
+        29137.947220361384,
+        48322.613544945292,  // ts 12 rejected
+        31056.413852819776,  // ts 12 repeated
+        37771.047066424137,
+        56236.288403836123,  // ts 14 rejected
+        54389.764270094922,  // ts 14 rejected
+        42756.662227525369,  // ts 14 repeated
+        60206.315291379709,
+        74747.111189489529,
+        114734.29990929153,
+        194708.67734889555,
+        354657.43222810357,
+        432000,
+        440640,
+        470880,
+        576720,
+        947160,  // ts 24 rejected
+        579960,  // ts 24 repeated
+        591300};
+
+    std::vector<double> const vec_t =
+        timeStepping(alg, nr_iterations, fixed_output_times, rejected_steps);
+
+    EXPECT_EQ(expected_vec_t.size(), vec_t.size());
+    ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                      5e4 * std::numeric_limits<double>::epsilon());
+}