Transport.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "Transport.h"

#include "algo/TransportEquations.h"
#include "math/OdeintFactory0.h"
#include "math/OdeintFactory2.h"
#include "math/OdeintTraits.h"

#include <boost/property_tree/ptree.hpp>
#include <cassert>

namespace SimPT_Default  {
namespace TimeEvolver {

using namespace std;
using namespace SimPT_Sim::ClockMan;

Transport::Transport(const CoreData& cd)
{
        Initialize(cd);
}

void Transport::Initialize(const CoreData& cd)
{
        assert( cd.Check() && "CoreData not ok.");
        m_cd = cd;
        const auto& p              = m_cd.m_parameters;

        m_abs_tolerance            = p->get<double>("ode_integration.abs_tolerance");
        m_is_stationary            = false;
        m_ode_solver               = p->get<string>("ode_integration.ode_solver", "runge_kutta_dopri5");
        m_rel_tolerance            = p->get<double>("ode_integration.rel_tolerance");
        m_small_time_increment     = p->get<double>("ode_integration.small_time_increment");
        m_stationarity_check       = p->get<bool>("termination.stationarity_check");
}

std::tuple<SimTimingTraits::CumulativeTimings, bool> Transport::operator()(double time_slice, SimPhase)
{
        // -------------------------------------------------------------------------------
        // Intro.
        // -------------------------------------------------------------------------------
        Stopclock chrono_stat("transport", true);
        bool is_stationary = false;
        TransportEquations equations(m_cd);

        // -------------------------------------------------------------------------------
        // Solve ODE's for cell2cell transport and diffusion.
        // -------------------------------------------------------------------------------
        const unsigned int num_eqns = equations.GetEquationCount();
        if (num_eqns > 0) {
                using State  = vector<double>;
                using System = OdeintTraits<State>::System;
                using Solver = OdeintTraits<State>::Solver;

                // Get current state of ODEs from the mesh
                State y_curr(num_eqns);
                equations.GetVariables(y_curr);

                // Wraps the ODE function dydt = F(y,t)
                System ode_function = [&equations](const State& y, State& dydt, const double) {
                        equations.SetVariables(y);
                        equations.GetDerivatives(dydt);
                };

                // Construct a suitable ODE solver either by looking at config or picking a default.
                Solver ode_solver;
                if (OdeintFactory0().ProvidesSolver(m_ode_solver)) {
                        // Solvers with fixed time steps
                        ode_solver = OdeintFactory0().Create(m_ode_solver);
                } else if (OdeintFactory2().ProvidesSolver(m_ode_solver)) {
                        // Solvers with adaptive time steps
                        ode_solver = OdeintFactory2().Create(m_ode_solver, m_abs_tolerance, m_rel_tolerance);
                } else {
                        // Pick a safe default in case none was given
                        ode_solver = OdeintFactory2().Create("runge_kutta_dopri5", m_abs_tolerance, m_rel_tolerance);
                }

                // ODE system's current start_time.
                const auto start_time = m_cd.m_time_data->m_sim_time;
                ode_solver(ode_function, y_curr, start_time, start_time + time_slice, m_small_time_increment);

                // Update mesh to reflect current state of ODEs
                equations.SetVariables(y_curr);

                // --------------------------------------------------------------------
                // Post-solve checks
                // --------------------------------------------------------------------
                if (m_stationarity_check) {
                        State dydt(num_eqns);
                        equations.GetDerivatives(dydt);
                        double norm = 0.0;
                        for (size_t i = 0; i < num_eqns; ++i) {
                                norm = max(norm, fabs(dydt[i]));
                        }
                        is_stationary = (norm <= m_abs_tolerance);
                }
        }

        // -------------------------------------------------------------------------------
        // Update simulated time.
        // -------------------------------------------------------------------------------
        m_cd.m_time_data->m_sim_time += time_slice;

        // -------------------------------------------------------------------------------
        // We are done ...
        // -------------------------------------------------------------------------------
        CumulativeTimings timings;
        timings.Record(chrono_stat.GetName(), chrono_stat.Get());
        return make_tuple(timings, is_stationary);
}

} // namespace
} // namespace