FixedTimeStepping.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "FixedTimeStepping.h"
 
#include <spdlog/fmt/ranges.h>
 
#include <algorithm>
#include <cassert>
#include <limits>
#include <numeric>
 
#include "NumLib/TimeStepping/Time.h"
 
namespace
{
NumLib::Time addTimeIncrement(std::vector<double>& delta_ts,
                              std::size_t const repeat, double const delta_t,
                              NumLib::Time const t_curr)
{
    auto const new_size = delta_ts.size() + repeat;
    try
    {
        delta_ts.resize(new_size, delta_t);
    }
    catch (std::exception const& e)
    {
        OGS_FATAL(
            "Resize of the time steps vector failed for the requested new size "
            "{:d}. Probably there is not enough memory ({:g} GiB "
            "requested).\nThrown exception: {:s}",
            new_size,
            new_size * sizeof(double) / 1024. / 1024. / 1024.,
            e.what());
    }
 
    // 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.
    return std::accumulate(end(delta_ts) - repeat, end(delta_ts), t_curr);
}
}  // namespace
 
namespace NumLib
{
std::size_t findDeltatInterval(Time const& t_initial,
                               std::vector<double> const& delta_ts,
                               Time 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];
    }
 
    return std::numeric_limits<std::size_t>::max();
}
 
void incorporateFixedTimesForOutput(
    NumLib::Time const t_initial, NumLib::Time 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,
            [](auto const time, NumLib::Time const& initial_time)
            { return NumLib::Time(time) < initial_time; });
        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, Time(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 = lower_bound + delta_ts[interval_number];
        // in order to use the comparison from struct Time
        if (NumLib::Time(fixed_time_for_output) == lower_bound)
        {
            continue;
        }
        if (upper_bound == NumLib::Time(fixed_time_for_output))
        {
            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);
    }
}
 
FixedTimeStepping::FixedTimeStepping(
    double t0, double tn, std::vector<RepeatDtPair> const& repeat_dt_pairs,
    std::vector<double> const& fixed_times_for_output)
    : TimeStepAlgorithm(t0, tn)
{
    Time t_curr = _t_initial;
 
    if (!areRepeatDtPairsValid(repeat_dt_pairs))
    {
        OGS_FATAL("FixedTimeStepping: Couldn't construct object from data");
    }
    for (auto const& [repeat, delta_t] : repeat_dt_pairs)
    {
        if (t_curr <= _t_end)
        {
            t_curr = addTimeIncrement(_dt_vector, repeat, delta_t, t_curr);
        }
    }
 
    // append last delta_t until t_end is reached
    if (t_curr <= _t_end)
    {
        auto const delta_t = std::get<1>(repeat_dt_pairs.back());
        auto const repeat = static_cast<std::size_t>(
            std::ceil((_t_end() - t_curr()) / delta_t));
        addTimeIncrement(_dt_vector, repeat, delta_t, t_curr);
    }
 
    incorporateFixedTimesForOutput(_t_initial, _t_end, _dt_vector,
                                   fixed_times_for_output);
}
 
FixedTimeStepping::FixedTimeStepping(double t0, double t_end, double dt)
    : TimeStepAlgorithm(t0, t_end)
{
    auto const new_size =
        static_cast<std::size_t>(std::ceil((t_end - t0) / dt));
    try
    {
        _dt_vector = std::vector<double>(new_size, dt);
    }
    catch (std::length_error const& e)
    {
        OGS_FATAL(
            "Resize of the time steps vector failed for the requested new "
            "size {}. Probably there is not enough memory ({:g} GiB "
            "requested).\n"
            "Thrown exception: {}",
            new_size, new_size * sizeof(double) / 1024. / 1024. / 1024.,
            e.what());
    }
    catch (std::bad_alloc const& e)
    {
        OGS_FATAL(
            "Allocation of the time steps vector failed for the requested "
            "size {}. Probably there is not enough memory ({:g} GiB "
            "requested).\n"
            "Thrown exception: {}",
            new_size,
            new_size * sizeof(double) / 1024. / 1024. / 1024.,
            e.what());
    }
}
 
double FixedTimeStepping::next(double const /*solution_error*/,
                               int const /*number_iterations*/,
                               NumLib::TimeStep& /*ts_previous*/,
                               NumLib::TimeStep& ts_current)
{
    // check if last time step
    if (ts_current.timeStepNumber() == _dt_vector.size() ||
        ts_current.current() >= end())
    {
        return 0.0;
    }
 
    double dt = _dt_vector[ts_current.timeStepNumber()];
    if (ts_current.current() + dt > end())
    {  // upper bound by t_end
        dt = end()() - ts_current.current()();
    }
 
    return dt;
}
 
bool FixedTimeStepping::areRepeatDtPairsValid(
    std::vector<RepeatDtPair> const& repeat_dt_pairs)
{
    if (repeat_dt_pairs.empty())
    {
        return false;
    }
 
    for (auto const& [repeat, delta_t] : repeat_dt_pairs)
    {
        if (repeat == 0)
        {
            ERR("FixedTimeStepping: <repeat> is zero.");
            return false;
        }
        if (delta_t <= 0.0)
        {
            ERR("FixedTimeStepping: timestep <delta_t> is <= 0.0.");
            return false;
        }
    }
    return true;
}
 
}  // namespace NumLib