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 0000000000000000000000000000000000000000..c112eb6756de21e4fb65a4ee7c71ac005fa336ed --- /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 506b4c4199d920ef0c54eb7d24658a9311fcd493..c57c891ceccacceaca8bc17406efdb9fd596d891 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/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp index 8529c65bd5a8af589c4daf7f6f217c09c0ab693f..1f8a692ac896a51a3b3cd64149557bf2345ef28c 100644 --- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp +++ b/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 diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h index 4819d27fa991dc513ef34fb5d4dea2a06c57fa68..af5a598e642cd72eb6fa23e7c7456a10fdf0af4f 100644 --- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h +++ b/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 diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp index 691b71b16470bb60b5164239924c634c03ce7ffe..a204a053b02c3621d184a1adb338d8d7c1e66a7b 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. 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; diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.h index 78ef8a2c097d3aa04e61084ccc1e9ded4b11167e..9dcf45eacbc2e783bd0ec650f69dab6cd246e144 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; @@ -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. diff --git a/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.cpp b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0f627b3bfc035153f41036d6279366a187e2cf60 --- /dev/null +++ b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.cpp @@ -0,0 +1,35 @@ +/** + * \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 diff --git a/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h index d0285855b379bf1410dc7448b6f3df3886b119cd..faf0a559f57379d1a7729ebd6300b5442d0b3b8f 100644 --- a/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h +++ b/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 diff --git a/Tests/Data/Parabolic/ComponentTransport/MassConservation/mass_conservation.prj b/Tests/Data/Parabolic/ComponentTransport/MassConservation/mass_conservation.prj index ac24f43a46088fd6ba7cd014f77c9d87ec56a32f..f24fa85dd8c641a4882ae361af6c71ad3cad94c3 100644 --- a/Tests/Data/Parabolic/ComponentTransport/MassConservation/mass_conservation.prj +++ b/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> diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp index 24a83dd91f8e58d1728cd75c26278c71e1dbc5f0..d9c9ab49fab1f7848762918b1a2a4ffa1878e84f 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()); +}