flow_simulation_base.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Rene Schmieding, Henrik Zunker
*
* 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 MIO_COMPARTMENTS_FLOW_SIMULATION_BASE_H
#define MIO_COMPARTMENTS_FLOW_SIMULATION_BASE_H
 
#include "memilio/compartments/flow_model.h"
#include "memilio/compartments/simulation_base.h"
#include "memilio/utils/time_series.h"
 
namespace mio
{
namespace details
{
 
template <typename FP, IsFlowModel<FP> M, class... Integrands>
class FlowSimulationBase : public SimulationBase<FP, M, Integrands...>
{
 
public:
    using Base  = SimulationBase<FP, M, Integrands...>;
    using Model = M;
    using Core  = IntegratorCore<FP, Integrands...>;
 
    FlowSimulationBase(Model const& model, std::unique_ptr<Core>&& integrator_core, FP t0, FP dt)
        : Base(model, std::move(integrator_core), t0, dt)
        , m_flow_result(t0, model.get_initial_flows())
    {
    }
 
    TimeSeries<FP>& get_flows()
    {
        return m_flow_result;
    }
 
    const TimeSeries<FP>& get_flows() const
    {
        return m_flow_result;
    }
protected:
    void compute_population_results()
    {
        const auto& flows = get_flows();
        const auto& model = this->get_model();
        auto& result      = this->get_result();
        // take the last time point as base result (instead of the initial results), so that we use external changes
        const size_t last_tp = result.get_num_time_points() - 1;
        // calculate new time points
        for (Eigen::Index i = result.get_num_time_points(); i < flows.get_num_time_points(); i++) {
            result.add_time_point(flows.get_time(i));
            model.get_derivatives(flows.get_value(i) - flows.get_value(last_tp), result.get_value(i));
            result.get_value(i) += result.get_value(last_tp);
        }
    }
 
private:
    mio::TimeSeries<FP> m_flow_result; 
};
 
template <typename FP, class M, class... Integrands>
class FlowSimulationBase<FP, M, SystemIntegrator<FP, Integrands...>> : public FlowSimulationBase<FP, M, Integrands...>
{
    using FlowSimulationBase<FP, M, Integrands...>::FlowSimulationBase;
};
 
} // namespace details
} // namespace mio
 
#endif // MIO_COMPARTMENTS_FLOW_SIMULATION_BASE_H