integrator.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef INTEGRATOR_H
#define INTEGRATOR_H
 
#include "memilio/config.h"
#include "memilio/math/floating_point.h"
#include "memilio/utils/time_series.h"
#include "memilio/utils/logging.h"
#include <memory>
#include <functional>
 
namespace mio
{
 
template <typename FP>
using DerivFunction =
    std::function<void(Eigen::Ref<const Eigen::VectorX<FP>> y, FP t, Eigen::Ref<Eigen::VectorX<FP>> dydt)>;
 
template <typename FP, class... Integrands>
class IntegratorCore
{
public:
    IntegratorCore(const FP& dt_min, const FP& dt_max)
        : m_dt_min(dt_min)
        , m_dt_max(dt_max)
    {
    }
 
    virtual ~IntegratorCore() {};
 
    virtual std::unique_ptr<IntegratorCore<FP, Integrands...>> clone() const = 0;
 
    virtual bool step(const Integrands&... fs, Eigen::Ref<const Eigen::VectorX<FP>> yt, FP& t, FP& dt,
                      Eigen::Ref<Eigen::VectorX<FP>> ytp1) const = 0;
 
    FP& get_dt_min()
    {
        return m_dt_min;
    }
    const FP& get_dt_min() const
    {
        return m_dt_min;
    }
    FP& get_dt_max()
    {
        return m_dt_max;
    }
    const FP& get_dt_max() const
    {
        return m_dt_max;
    }
private:
    FP m_dt_min, m_dt_max; 
};
 
template <class FP>
using OdeIntegratorCore = IntegratorCore<FP, DerivFunction<FP>>;
 
template <class FP>
using SdeIntegratorCore = IntegratorCore<FP, DerivFunction<FP>, DerivFunction<FP>>;
 
template <typename FP, class... Integrands>
class SystemIntegrator
{
public:
    using Core = IntegratorCore<FP, Integrands...>;
    SystemIntegrator(std::unique_ptr<Core>&& core)
        : m_core(std::move(core))
        , m_is_adaptive(false)
    {
    }
 
    SystemIntegrator(const SystemIntegrator& other)
        : m_core(other.m_core->clone())
        , m_is_adaptive(other.m_is_adaptive)
    {
    }
 
    SystemIntegrator& operator=(const SystemIntegrator& other)
    {
        if (this != &other) {
            m_core        = other.m_core->clone();
            m_is_adaptive = other.m_is_adaptive;
        }
        return *this;
    }
 
    Eigen::Ref<Eigen::VectorX<FP>> advance(const Integrands&... fs, const FP tmax, FP& dt, TimeSeries<FP>& results)
    {
        // hint at std functions for ADL
        using std::ceil;
        using std::fabs;
        using std::max;
        using std::min;
        const FP t0  = results.get_last_time();
        const FP eps = Limits<FP>::zero_tolerance();
        assert(tmax > t0);
        assert(dt > 0);
 
        const auto num_steps =
            static_cast<size_t>(ceil((tmax - t0) / dt)); // estimated number of time steps (if equidistant)
 
        results.reserve(results.get_num_time_points() + num_steps);
 
        bool step_okay = true;
 
        FP dt_copy; // used to check whether step sizing is adaptive
        FP dt_restore     = 0.0; // used to restore dt, if dt was decreased to reach tmax
        FP dt_min_restore = m_core->get_dt_min(); // used to restore dt_min, if it was decreased to reach tmax
        FP t              = t0;
 
        for (size_t i = results.get_num_time_points() - 1; fabs((tmax - t) / (tmax - t0)) > eps; ++i) {
            // We don't make time steps too small as the error estimator of an adaptive integrator
            //may not be able to handle it. this is very conservative and maybe unnecessary,
            //but also unlikely to happen. may need to be reevaluated.
 
            if (dt > tmax - t) {
                dt_restore = dt;
                dt         = tmax - t;
                // if necessary, also reduce minimal step size such that we do not step past tmax
                m_core->get_dt_min() = min<FP>(tmax - t, m_core->get_dt_min());
                // if dt_min was reduced, the following step will be the last due to dt == dt_min (see step method)
                // dt_min must be restored after this loop
            }
 
            dt_copy = dt;
 
            results.add_time_point();
            step_okay &= m_core->step(fs..., results[i], t, dt, results[i + 1]);
            results.get_last_time() = t;
 
            // if dt has been changed by step, register the current m_core as adaptive.
            m_is_adaptive |= !floating_point_equal<FP>(dt, dt_copy, eps);
        }
        m_core->get_dt_min() = dt_min_restore; // restore dt_min
        // if dt was decreased to reach tmax in the last time iteration,
        // we restore it as it is now probably smaller than required for tolerances
        dt = max<FP>(dt, dt_restore);
 
        if (m_is_adaptive) {
            if (!step_okay) {
                log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
            }
            else if (fabs((tmax - t) / (tmax - t0)) > eps) {
                log_warning("Last time step too small. Could not reach tmax exactly.");
            }
            else {
                log_info("Adaptive step sizing successful to tolerances.");
            }
        }
 
        return results.get_last_value();
    }
 
    void set_integrator_core(std::unique_ptr<Core>&& core)
    {
        m_core.swap(core);
        m_is_adaptive = false;
    }
 
    Core& get_integrator_core()
    {
        return *m_core;
    }
    const Core& get_integrator_core() const
    {
        return *m_core;
    }
 
private:
    std::unique_ptr<Core> m_core;
    bool m_is_adaptive;
};
 
template <typename FP>
using OdeIntegrator = SystemIntegrator<FP, DerivFunction<FP>>;
 
template <typename FP>
using SdeIntegrator = SystemIntegrator<FP, DerivFunction<FP>, DerivFunction<FP>>;
 
} // namespace mio
 
#endif // INTEGRATOR_H