From d38bb0ca80eca04a2a3af20422702973e08f77ab Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 1 Dec 2022 15:10:50 +0100
Subject: [PATCH] [NL] Include fixed output times in FixedTimeStepping
algorithm
---
.../Algorithms/CreateFixedTimeStepping.cpp | 119 ++++++++++++++++--
.../Algorithms/CreateFixedTimeStepping.h | 4 +
2 files changed, 115 insertions(+), 8 deletions(-)
diff --git a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
index 9d6c8bec076..461ec71f28d 100644
--- a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
@@ -11,8 +11,13 @@
#include "CreateFixedTimeStepping.h"
+#include <fmt/ranges.h>
+
+#include <algorithm>
+#include <numeric>
#include <string>
+#include "BaseLib/Algorithm.h"
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Error.h"
#include "FixedTimeStepping.h"
@@ -20,6 +25,102 @@
namespace NumLib
{
+std::size_t findDeltatInterval(double const t_initial,
+ std::vector<double> const& delta_ts,
+ double const fixed_output_time)
+{
+ if (fixed_output_time < t_initial)
+ {
+ return std::numeric_limits<std::size_t>::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<std::size_t>::max();
+ }
+ return std::numeric_limits<std::size_t>::max();
+}
+
+void incorporateFixedTimesForOutput(
+ double const t_initial, double const t_end, std::vector<double>& delta_ts,
+ std::vector<double> 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<std::size_t>::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<double>::epsilon())
+ {
+ continue;
+ }
+ if (upper_bound - fixed_time_for_output <=
+ std::numeric_limits<double>::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<TimeStepAlgorithm> createFixedTimeStepping(
BaseLib::ConfigTree const& config,
@@ -33,15 +134,15 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
//! \ogs_file_param{prj__time_loop__processes__process__time_stepping__FixedTimeStepping__t_end}
auto const t_end = config.getConfigParameter<double>("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<double> timesteps;
+ std::vector<double> 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");
@@ -64,10 +165,10 @@ std::unique_ptr<TimeStepAlgorithm> 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)
{
@@ -100,10 +201,10 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
{
auto const repeat =
static_cast<std::size_t>(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)
{
@@ -128,6 +229,8 @@ std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
}
}
- return std::make_unique<FixedTimeStepping>(t_initial, t_end, timesteps);
+ incorporateFixedTimesForOutput(t_initial, t_end, delta_ts,
+ fixed_times_for_output);
+ return std::make_unique<FixedTimeStepping>(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 fd0aa7fe8d2..9c2062faf86 100644
--- a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h
+++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h
@@ -23,6 +23,10 @@ namespace NumLib
{
class TimeStepAlgorithm;
+std::size_t findDeltatInterval(double const t_initial,
+ std::vector<double> const& delta_ts,
+ double const fixed_output_time);
+
/// Create a FixedTimeStepping time stepper from the given
/// configuration
std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
--
GitLab