diff --git a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp index b80498fee37..c88d5862eef 100644 --- a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp +++ b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp @@ -20,7 +20,8 @@ namespace NumLib { class TimeStepAlgorithm; std::unique_ptr createEvolutionaryPIDcontroller( - BaseLib::ConfigTree const& config) + BaseLib::ConfigTree const& config, + std::vector 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 createEvolutionaryPIDcontroller( //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__EvolutionaryPIDcontroller__tol} auto const tol = config.getConfigParameter("tol"); - return std::make_unique( - t0, t_end, h0, h_min, h_max, rel_h_min, rel_h_max, tol); + return std::make_unique(t0, + t_end, + h0, + h_min, + h_max, + rel_h_min, + rel_h_max, + tol, + fixed_times_for_output); } } // end of namespace NumLib diff --git a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h index f22181680ac..3a38a86a5ca 100644 --- a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h +++ b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h @@ -12,6 +12,7 @@ #pragma once #include +#include namespace BaseLib { @@ -25,5 +26,6 @@ class TimeStepAlgorithm; /// Create an EvolutionaryPIDcontroller time stepper from the given /// configuration std::unique_ptr createEvolutionaryPIDcontroller( - BaseLib::ConfigTree const& config); + BaseLib::ConfigTree const& config, + std::vector const& fixed_times_for_output); } // end of namespace NumLib diff --git a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp index ff7c7c6d19a..82e56c6103e 100644 --- a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp +++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp @@ -11,8 +11,13 @@ #include "CreateFixedTimeStepping.h" +#include + +#include +#include #include +#include "BaseLib/Algorithm.h" #include "BaseLib/ConfigTree.h" #include "BaseLib/Error.h" #include "FixedTimeStepping.h" @@ -20,9 +25,106 @@ namespace NumLib { +std::size_t findDeltatInterval(double const t_initial, + std::vector const& delta_ts, + double const fixed_output_time) +{ + if (fixed_output_time < t_initial) + { + return std::numeric_limits::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::max(); + } + return std::numeric_limits::max(); +} + +void incorporateFixedTimesForOutput( + double const t_initial, double const t_end, std::vector& delta_ts, + std::vector 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::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::epsilon()) + { + continue; + } + if (upper_bound - fixed_time_for_output <= + std::numeric_limits::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 createFixedTimeStepping( - BaseLib::ConfigTree const& config) + BaseLib::ConfigTree const& config, + std::vector 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 createFixedTimeStepping( //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__FixedTimeStepping__t_end} auto const t_end = config.getConfigParameter("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 timesteps; + std::vector 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 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 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 createFixedTimeStepping( { auto const repeat = static_cast(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 createFixedTimeStepping( } } - return std::make_unique(t_initial, t_end, timesteps); + incorporateFixedTimesForOutput(t_initial, t_end, delta_ts, + fixed_times_for_output); + return std::make_unique(t_initial, t_end, delta_ts); } } // end of namespace NumLib diff --git a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h index 444888b0fe6..9c2062faf86 100644 --- a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h +++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h @@ -12,6 +12,7 @@ #pragma once #include +#include namespace BaseLib { @@ -22,8 +23,13 @@ namespace NumLib { class TimeStepAlgorithm; +std::size_t findDeltatInterval(double const t_initial, + std::vector const& delta_ts, + double const fixed_output_time); + /// Create a FixedTimeStepping time stepper from the given /// configuration std::unique_ptr createFixedTimeStepping( - BaseLib::ConfigTree const& config); + BaseLib::ConfigTree const& config, + std::vector const& fixed_times_for_output); } // end of namespace NumLib diff --git a/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp index 8373e967fb1..3e3467280b6 100644 --- a/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp +++ b/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.cpp @@ -20,7 +20,8 @@ namespace NumLib { class TimeStepAlgorithm; std::unique_ptr createIterationNumberBasedTimeStepping( - BaseLib::ConfigTree const& config) + BaseLib::ConfigTree const& config, + std::vector 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 createIterationNumberBasedTimeStepping( return std::make_unique( 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 diff --git a/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.h b/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.h index a5a224a08bc..db45c66b156 100644 --- a/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.h +++ b/NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.h @@ -10,6 +10,7 @@ #pragma once #include +#include namespace BaseLib { @@ -26,5 +27,6 @@ namespace NumLib /// Create a IterationNumberBasedTimeStepping time stepper from the given /// configuration. std::unique_ptr createIterationNumberBasedTimeStepping( - BaseLib::ConfigTree const& config); + BaseLib::ConfigTree const& config, + std::vector const& fixed_times_for_output); } // namespace NumLib diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp index 6c37736e509..2856921ec6c 100644 --- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp +++ b/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), + [×tep_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::epsilon() * + timestep_current.current()) + { + return *fixed_output_time_it - timestep_current.current(); + } + } + } return limited_h; } diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h index 9238b3e4e86..c267eb442fd 100644 --- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h +++ b/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 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 const _fixed_times_for_output; + /** * Force the computed time step size in the given range * (see the formulas in the documentation of the class) diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp index ca1115bba7b..ac2ca65fca8 100644 --- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp +++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp @@ -26,14 +26,16 @@ IterationNumberBasedTimeStepping::IterationNumberBasedTimeStepping( double const t_initial, double const t_end, double const min_dt, double const max_dt, double const initial_dt, std::vector&& iter_times_vector, - std::vector&& multiplier_vector) + std::vector&& multiplier_vector, + std::vector const& fixed_times_for_output) : TimeStepAlgorithm(t_initial, t_end), _iter_times_vector(std::move(iter_times_vector)), _multiplier_vector(std::move(multiplier_vector)), _min_dt(min_dt), _max_dt(max_dt), _initial_dt(initial_dt), - _max_iter(_iter_times_vector.empty() ? 0 : _iter_times_vector.back()) + _max_iter(_iter_times_vector.empty() ? 0 : _iter_times_vector.back()), + _fixed_times_for_output(fixed_times_for_output) { if (_iter_times_vector.empty()) { @@ -135,6 +137,32 @@ double IterationNumberBasedTimeStepping::getNextTimeStepSize( dt = ts_previous.dt() * findMultiplier(_iter_times, ts_current); } + if (_fixed_times_for_output.empty()) + { + return std::clamp(dt, _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( + std::begin(_fixed_times_for_output), std::end(_fixed_times_for_output), + [&ts_current](auto const fixed_output_time) + { return ts_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 + dt) + if (*fixed_output_time_it < ts_current.current() + dt) + { + // check if the potential adjusted time step is larger than zero + if (std::abs(*fixed_output_time_it - ts_current.current()) > + std::numeric_limits::epsilon() * ts_current.current()) + { + return *fixed_output_time_it - ts_current.current(); + } + } + } return std::clamp(dt, _min_dt, _max_dt); } diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h index a6952f0cc81..5bfd3d65b9a 100644 --- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h +++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h @@ -81,14 +81,17 @@ class IterationNumberBasedTimeStepping final : public TimeStepAlgorithm * (\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_times_for_output a vector of fixed time points for output */ - IterationNumberBasedTimeStepping(double const t_initial, - double const t_end, - double const min_dt, - double const max_dt, - double const initial_dt, - std::vector&& iter_times_vector, - std::vector&& multiplier_vector); + IterationNumberBasedTimeStepping( + double const t_initial, + double const t_end, + double const min_dt, + double const max_dt, + double const initial_dt, + std::vector&& iter_times_vector, + std::vector&& multiplier_vector, + std::vector const& fixed_times_for_output); ~IterationNumberBasedTimeStepping() override = default; @@ -128,6 +131,7 @@ class IterationNumberBasedTimeStepping final : public TimeStepAlgorithm int _iter_times = 0; bool _previous_time_step_accepted = true; + std::vector const _fixed_times_for_output; }; } // namespace NumLib diff --git a/NumLib/TimeStepping/CreateTimeStepper.cpp b/NumLib/TimeStepping/CreateTimeStepper.cpp index 993582ef8ff..3da7bd9ad17 100644 --- a/NumLib/TimeStepping/CreateTimeStepper.cpp +++ b/NumLib/TimeStepping/CreateTimeStepper.cpp @@ -24,7 +24,8 @@ namespace NumLib { std::unique_ptr createTimeStepper( - BaseLib::ConfigTree const& config) + BaseLib::ConfigTree const& config, + std::vector const& fixed_times_for_output) { //! \ogs_file_param{prj__time_loop__processes__process__time_stepping__type} auto const type = config.peekConfigParameter("type"); @@ -37,15 +38,17 @@ std::unique_ptr createTimeStepper( } if (type == "FixedTimeStepping") { - return NumLib::createFixedTimeStepping(config); + return NumLib::createFixedTimeStepping(config, fixed_times_for_output); } if (type == "EvolutionaryPIDcontroller") { - return NumLib::createEvolutionaryPIDcontroller(config); + return NumLib::createEvolutionaryPIDcontroller(config, + fixed_times_for_output); } if (type == "IterationNumberBasedTimeStepping") { - return NumLib::createIterationNumberBasedTimeStepping(config); + return NumLib::createIterationNumberBasedTimeStepping( + config, fixed_times_for_output); } OGS_FATAL( "Unknown time stepping type: '{:s}'. The available types are: " diff --git a/NumLib/TimeStepping/CreateTimeStepper.h b/NumLib/TimeStepping/CreateTimeStepper.h index bf582927bac..2b3d45ee61f 100644 --- a/NumLib/TimeStepping/CreateTimeStepper.h +++ b/NumLib/TimeStepping/CreateTimeStepper.h @@ -12,6 +12,7 @@ #pragma once #include +#include namespace BaseLib { @@ -22,5 +23,6 @@ namespace NumLib { class TimeStepAlgorithm; std::unique_ptr createTimeStepper( - BaseLib::ConfigTree const& config); + BaseLib::ConfigTree const& config, + std::vector const& fixed_times_for_output); }; // namespace NumLib diff --git a/ProcessLib/CreateProcessData.cpp b/ProcessLib/CreateProcessData.cpp index f0af74ac6b1..cde92a31db7 100644 --- a/ProcessLib/CreateProcessData.cpp +++ b/ProcessLib/CreateProcessData.cpp @@ -74,7 +74,8 @@ std::vector> createPerProcessData( std::vector> const& processes, std::map> const& nonlinear_solvers, - bool const compensate_non_equilibrium_initial_residuum) + bool const compensate_non_equilibrium_initial_residuum, + std::vector const& fixed_times_for_output) { std::vector> per_process_data; int process_id = 0; @@ -103,7 +104,8 @@ std::vector> createPerProcessData( auto timestepper = NumLib::createTimeStepper( //! \ogs_file_param{prj__time_loop__processes__process__time_stepping} - pcs_config.getConfigSubtree("time_stepping")); + pcs_config.getConfigSubtree("time_stepping"), + fixed_times_for_output); auto conv_crit = NumLib::createConvergenceCriterion( //! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion} diff --git a/ProcessLib/CreateProcessData.h b/ProcessLib/CreateProcessData.h index f00cd5b1535..1c489b0dca0 100644 --- a/ProcessLib/CreateProcessData.h +++ b/ProcessLib/CreateProcessData.h @@ -19,6 +19,7 @@ std::vector> createPerProcessData( std::vector> const& processes, std::map> const& nonlinear_solvers, - bool const compensate_non_equilibrium_initial_residuum); + bool const compensate_non_equilibrium_initial_residuum, + std::vector const& fixed_times_for_output); } // namespace ProcessLib diff --git a/ProcessLib/CreateTimeLoop.cpp b/ProcessLib/CreateTimeLoop.cpp index c6727f4d142..9b1ada6225a 100644 --- a/ProcessLib/CreateTimeLoop.cpp +++ b/ProcessLib/CreateTimeLoop.cpp @@ -69,6 +69,8 @@ std::unique_ptr createTimeLoop( //! \ogs_file_param{prj__time_loop__outputs} : createOutputs(config.getConfigSubtree("outputs"), output_directory, meshes); + auto const fixed_times_for_output = + calculateUniqueFixedTimesForAllOutputs(outputs); if (auto const submesh_residuum_output_config_tree = //! \ogs_file_param{prj__time_loop__submesh_residuum_output} @@ -104,7 +106,7 @@ std::unique_ptr createTimeLoop( auto per_process_data = createPerProcessData( //! \ogs_file_param{prj__time_loop__processes} config.getConfigSubtree("processes"), processes, nonlinear_solvers, - compensate_non_equilibrium_initial_residuum); + compensate_non_equilibrium_initial_residuum, fixed_times_for_output); const bool use_staggered_scheme = ranges::any_of(processes.begin(), processes.end(), diff --git a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp index deac97076c9..db9ef994f1c 100644 --- a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp +++ b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp @@ -28,7 +28,7 @@ std::unique_ptr createTestTimeStepper( BaseLib::ConfigTree conf(std::move(ptree), "", BaseLib::ConfigTree::onerror, BaseLib::ConfigTree::onwarning); auto const& sub_config = conf.getConfigSubtree("time_stepping"); - return NumLib::createEvolutionaryPIDcontroller(sub_config); + return NumLib::createEvolutionaryPIDcontroller(sub_config, {}); } TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller) diff --git a/Tests/NumLib/TestTimeSteppingFixed.cpp b/Tests/NumLib/TestTimeSteppingFixed.cpp index 18d67fc806a..e3ba1b581c0 100644 --- a/Tests/NumLib/TestTimeSteppingFixed.cpp +++ b/Tests/NumLib/TestTimeSteppingFixed.cpp @@ -12,68 +12,320 @@ #include +#include #include +#include "NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h" #include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h" #include "NumLib/TimeStepping/TimeStep.h" #include "Tests/TestTools.h" #include "TimeSteppingTestingTools.h" -TEST(NumLib, TimeSteppingFixed) +namespace NumLib { - std::vector const dummy_number_iterations = {0, 0, 0, 0}; - // homogeneous dt - { - NumLib::FixedTimeStepping fixed(1, 31, 10); - const std::vector expected_vec_t = {1, 11, 21, 31}; +extern void incorporateFixedTimesForOutput( + double t_initial, double t_end, std::vector& timesteps, + std::vector const& fixed_times_for_output); +} - std::vector vec_t = - timeStepping(fixed, dummy_number_iterations, {}); +class NumLibTimeSteppingFixed_FindDeltaT : public ::testing::Test +{ + std::vector dts; - ASSERT_EQ(expected_vec_t.size(), vec_t.size()); - ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), - std::numeric_limits::epsilon()); +protected: + double const initial_dt = 1e-12; + + void SetUp() override + { + dts = {initial_dt}; + dts.insert(dts.end(), 10, 1e-1); } - // dt vector (t_end == t0 + sum(dt)) + std::size_t find(double const t) { - const std::vector fixed_dt = {10, 10, 10}; - NumLib::FixedTimeStepping fixed(1, 31, fixed_dt); - const std::vector expected_vec_t = {1, 11, 21, 31}; + double const t_initial = 0.0; + return NumLib::findDeltatInterval(t_initial, dts, t); + } +}; - std::vector vec_t = - timeStepping(fixed, dummy_number_iterations, {}); +TEST_F(NumLibTimeSteppingFixed_FindDeltaT, TestPointWithinMiddleInterval) +{ + EXPECT_EQ(5, find(0.5)); +} - ASSERT_EQ(expected_vec_t.size(), vec_t.size()); - ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), - std::numeric_limits::epsilon()); - } +TEST_F(NumLibTimeSteppingFixed_FindDeltaT, TestPointAtBoundary) +{ + EXPECT_EQ(5, find(0.5 + std::numeric_limits::epsilon())); +} - // dt vector (t_end < t0 + sum(dt)) - { - const std::vector fixed_dt = {5, 10, 20}; - NumLib::FixedTimeStepping fixed(1, 31, fixed_dt); - const std::vector expected_vec_t = {1, 6, 16, 31}; +TEST_F(NumLibTimeSteppingFixed_FindDeltaT, TestPointJustBeforeNextInterval) +{ + EXPECT_EQ(6, find(0.6 - initial_dt)); +} + +TEST_F(NumLibTimeSteppingFixed_FindDeltaT, TestPointBeforeInitialTime) +{ + EXPECT_EQ(std::numeric_limits::max(), find(-initial_dt)); +} + +TEST_F(NumLibTimeSteppingFixed_FindDeltaT, TestPointBeyondEndTime) +{ + EXPECT_EQ(std::numeric_limits::max(), find(1000)); +} + +class NumLibTimeSteppingFixed_TimeSteps : public ::testing::Test +{ + std::vector const dummy_number_iterations = {0, 0, 0, 0}; - std::vector vec_t = - timeStepping(fixed, dummy_number_iterations, {}); +protected: + void checkExpectedTimes(NumLib::FixedTimeStepping& algorithm, + std::vector const expected_times) const + { + std::vector times = + timeStepping(algorithm, dummy_number_iterations, {}); - ASSERT_EQ(expected_vec_t.size(), vec_t.size()); - ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), + ASSERT_EQ(expected_times.size(), times.size()); + ASSERT_ARRAY_NEAR(expected_times, times, expected_times.size(), std::numeric_limits::epsilon()); } +}; - // dt vector (t_end > t0 + sum(dt)) - { - const std::vector fixed_dt = {5, 10, 10}; - NumLib::FixedTimeStepping fixed(1, 31, fixed_dt); - const std::vector expected_vec_t = {1, 6, 16, 26}; +TEST_F(NumLibTimeSteppingFixed_TimeSteps, HomogeneousDt) +{ + NumLib::FixedTimeStepping algorithm(1, 31, 10); + + checkExpectedTimes(algorithm, {1, 11, 21, 31}); +} - std::vector vec_t = - timeStepping(fixed, dummy_number_iterations, {}); +TEST_F(NumLibTimeSteppingFixed_TimeSteps, DtVectorEqualSum) +{ + std::vector const fixed_dt = {10, 10, 10}; + NumLib::FixedTimeStepping algorithm(1, 31, fixed_dt); - ASSERT_EQ(expected_vec_t.size(), vec_t.size()); - ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), - std::numeric_limits::epsilon()); + checkExpectedTimes(algorithm, {1, 11, 21, 31}); +} + +TEST_F(NumLibTimeSteppingFixed_TimeSteps, DtVectorLessThanSum) +{ + std::vector const fixed_dt = {5, 10, 20}; + NumLib::FixedTimeStepping algorithm(1, 31, fixed_dt); + + checkExpectedTimes(algorithm, {1, 6, 16, 31}); +} + +TEST_F(NumLibTimeSteppingFixed_TimeSteps, DtVectorGreaterThanSum) +{ + std::vector const fixed_dt = {5, 10, 10}; + NumLib::FixedTimeStepping algorithm(1, 31, fixed_dt); + + checkExpectedTimes(algorithm, {1, 6, 16, 26}); +} + +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, incorporateFixedTimesForOutput_2) +{ + double t_initial = 1e-10; + std::vector dts{}; + dts.insert(dts.end(), 10, 1e-1); + std::vector fixed_times_for_output{{0.5, 1}}; + + auto const expected_time = std::accumulate(dts.begin(), dts.end(), 0.0); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, dts, + fixed_times_for_output); + + // incorporation of time steps doesn't influence the entire simulation time + EXPECT_NEAR(expected_time, std::accumulate(dts.begin(), dts.end(), 0.0), + std::numeric_limits::epsilon()); + + ASSERT_EQ(1e-1, dts[0]); + ASSERT_EQ(1e-1, dts[1]); + ASSERT_EQ(1e-1, dts[2]); + ASSERT_EQ(1e-1, dts[3]); // time point 0.4 + 1e-10 + EXPECT_NEAR(1e-1 - t_initial, dts[4], + std::numeric_limits::epsilon()); // time point 0.5 + EXPECT_EQ(5e-1 + 1e-10 - 0.5, dts[5]); // time point 0.5 + 1e-10 + EXPECT_EQ(1e-1, dts[6]); // time point 0.6 + 1e-10 + EXPECT_EQ(1e-1, dts[7]); // time point 0.7 + 1e-10 + EXPECT_EQ(1e-1, dts[8]); // time point 0.8 + 1e-10 + EXPECT_EQ(1e-1, dts[9]); // time point 0.9 + 1e-10 + EXPECT_NEAR(1.0 - (9e-1 + t_initial), dts[10], + std::numeric_limits::epsilon()); // time point 1.0 + EXPECT_NEAR( + t_initial, dts[11], + std::numeric_limits::epsilon()); // time point 1.0 + 1e-10 +} + +// unit test related to ThermoRichardsMechanics/LiakopoulosHM/liakopoulos.prj +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, incorporateFixedTimesForOutput_3) +{ + double t_initial = 0.0; + std::vector dts{}; + + // 10 1 + dts.insert(dts.end(), 10, 1); + // 9 10 + dts.insert(dts.end(), 9, 10); + // 11 100 + dts.insert(dts.end(), 11, 100); + // 3 200 + dts.insert(dts.end(), 3, 200); + // 3 400 + dts.insert(dts.end(), 3, 400); + // 7 600 + dts.insert(dts.end(), 7, 600); + + std::vector fixed_times_for_output{{0.06, 60., 120., 300.0, 600.0, + 1200.0, 2400.0, 4800.0, 6000.0, + 7200.0}}; + + auto const expected_time = std::accumulate(dts.begin(), dts.end(), 0.0); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, dts, + fixed_times_for_output); + + // incorporation of time steps doesn't influence the entire simulation time + ASSERT_EQ(expected_time, std::accumulate(dts.begin(), dts.end(), 0.0)); + + ASSERT_EQ(0.06, dts[0]); + ASSERT_EQ(1 - 0.06, dts[1]); + for (std::size_t k = 2; k < 11; ++k) + { + ASSERT_EQ(1, dts[k]); + } + for (std::size_t k = 11; k < 20; ++k) + { + ASSERT_EQ(10, dts[k]); } + ASSERT_EQ(20, dts[20]); // step size of 100 is split into 20 + 80 because + // of output at 120 + ASSERT_EQ(80, dts[21]); + for (std::size_t k = 22; k < 32; ++k) + { + ASSERT_EQ(100, dts[k]); + } + for (std::size_t k = 32; k < 35; ++k) + { + ASSERT_EQ(200, dts[k]); + } + ASSERT_EQ(400, dts[35]); + ASSERT_EQ(200, dts[36]); // step size of 400 is split into 200 + 200 + // because of output at 2400 + ASSERT_EQ(200, dts[37]); + ASSERT_EQ(400, dts[38]); + for (std::size_t k = 39; k < 46; ++k) + { + ASSERT_EQ(600, dts[k]); + } +} + +// unit test related to +// ogs-ThermoMechanics_CreepBGRa_Verification_m2_3Dload_m2_3Dload +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, incorporateFixedTimesForOutput_4) +{ + double t_initial = 0.0; + std::vector dts{t_initial}; + + // 1 1e-10 + dts.insert(dts.end(), 1, 1e-10); + // 4 0.01 + dts.insert(dts.end(), 4, 1e-2); + // 1 0.0099999999 + dts.insert(dts.end(), 1, 0.0099999999); + // 100 0.01 + dts.insert(dts.end(), 100, 1e-2); + + std::vector fixed_times_for_output{{0.5, 1.}}; + + auto const expected_end_time = + std::accumulate(dts.begin(), dts.end(), t_initial); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_end_time, dts, + fixed_times_for_output); + + // incorporation of time steps doesn't influence the entire simulation time + ASSERT_EQ(expected_end_time, + std::accumulate(dts.begin(), dts.end(), t_initial)); +} + +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, incorporateFixedTimesForOutput) +{ + double t_initial = 1.0; + std::vector dts{{10, 10, 10}}; + std::vector fixed_times_for_output{{9, 12, 28}}; + + auto const expected_time = std::accumulate(dts.begin(), dts.end(), 0.0); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, dts, + fixed_times_for_output); + + ASSERT_EQ(expected_time, std::accumulate(dts.begin(), dts.end(), 0.0)); + ASSERT_EQ(8.0, dts[0]); + ASSERT_EQ(2.0, dts[1]); + ASSERT_EQ(1.0, dts[2]); + ASSERT_EQ(9.0, dts[3]); + ASSERT_EQ(7.0, dts[4]); + ASSERT_EQ(3.0, dts[5]); +} + +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, + incorporateFixedTimesForOutput_Matching) +{ + double t_initial = 1.0; + std::vector timesteps{{10, 10, 10}}; + std::vector fixed_times_for_output{{9, 11, 31}}; + + auto const expected_time = + std::accumulate(timesteps.begin(), timesteps.end(), 0.0); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, timesteps, + fixed_times_for_output); + + ASSERT_EQ(expected_time, + std::accumulate(timesteps.begin(), timesteps.end(), 0.0)); + ASSERT_EQ(8.0, timesteps[0]); + ASSERT_EQ(2.0, timesteps[1]); + ASSERT_EQ(10.0, timesteps[2]); + ASSERT_EQ(10.0, timesteps[3]); +} + +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, + OutputTimeBeforeSimulationStartTime) +{ + double t_initial = 10.0; + std::vector timesteps{{10, 10, 10}}; + std::vector fixed_times_for_output{{9, 12, 28}}; + + auto const expected_time = + std::accumulate(timesteps.begin(), timesteps.end(), t_initial); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, timesteps, + fixed_times_for_output); + ASSERT_EQ(expected_time, + std::accumulate(timesteps.begin(), timesteps.end(), t_initial)); + ASSERT_EQ(2.0, timesteps[0]); + ASSERT_EQ(8.0, timesteps[1]); + ASSERT_EQ(8.0, timesteps[2]); + ASSERT_EQ(2.0, timesteps[3]); + ASSERT_EQ(10.0, timesteps[4]); +} + +TEST(NumLibTimeSteppingFixed_FixedOutputTimes, OutputTimeAfterSimulationEndTime) +{ + double t_initial = 1.0; + std::vector timesteps{{10, 10, 10}}; + std::vector fixed_times_for_output{{9, 12, 28, 33}}; + + auto const expected_time = + std::accumulate(timesteps.begin(), timesteps.end(), t_initial); + + NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, timesteps, + fixed_times_for_output); + ASSERT_EQ(expected_time, + std::accumulate(timesteps.begin(), timesteps.end(), t_initial)); + ASSERT_EQ(8.0, timesteps[0]); + ASSERT_EQ(2.0, timesteps[1]); + ASSERT_EQ(1.0, timesteps[2]); + ASSERT_EQ(9.0, timesteps[3]); + ASSERT_EQ(7.0, timesteps[4]); + ASSERT_EQ(3.0, timesteps[5]); } diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp index e7681be237a..5de45111281 100644 --- a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp +++ b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp @@ -25,9 +25,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased1) { std::vector iter_times_vector = {0, 3, 5, 7}; std::vector 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.; const double end_time = alg.end(); @@ -145,9 +145,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased2) { std::vector iter_times_vector = {0, 3, 5, 7}; std::vector 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 nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1}; const std::vector expected_vec_t = {1, 2, 4, 8, 16, @@ -165,9 +165,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased2FixedOutputTimes) std::vector iter_times_vector = {0, 3, 5, 7}; std::vector multiplier_vector = {2.0, 1.0, 0.5, 0.25}; std::vector fixed_output_times = {5, 20}; - 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 nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1, 1, 1, 1, 1}; const std::vector expected_vec_t = {1, 2, 4, 5, 7, 9, 10,