limex_integrator.hpp

/*
 * SPDX-FileCopyrightText: Copyright (c) 2014-2025:  G-CSC, Goethe University Frankfurt
 * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
 *
 * Author: Arne Naegel
 *
 * This file is part of UG4.
 *
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 *
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 *
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 *
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef LIMEX_INTEGRATOR_HPP_
#define LIMEX_INTEGRATOR_HPP_
/*
#define XMTHREAD_BOOST
#ifdef XMTHREAD_BOOST
#include <boost/thread/thread.hpp>
#endif
*/
#include <string>
 
#include "common/stopwatch.h"
 
#include "lib_algebra/operator/interface/operator.h"
#include "lib_algebra/operator/interface/operator_inverse.h"
#include "lib_algebra/operator/linear_solver/linear_solver.h"
#include "lib_algebra/operator/debug_writer.h"
 
#include "lib_disc/function_spaces/grid_function.h"
#include "lib_disc/assemble_interface.h" // TODO: missing IAssemble in following file:
#include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
#include "lib_disc/operator/non_linear_operator/assembled_non_linear_operator.h"
#include "lib_disc/spatial_disc/domain_disc.h"
#include "lib_disc/time_disc/time_disc_interface.h"
#include "lib_disc/time_disc/theta_time_step.h"
#include "lib_disc/time_disc/solution_time_series.h"
#include "lib_disc/function_spaces/grid_function_util.h" // SaveVectorForConnectionViewer
#include "lib_disc/function_spaces/metric_spaces.h"
#include "lib_disc/io/vtkoutput.h"
 
#include "lib_grid/refinement/refiner_interface.h"
 
 
// own headers
#include "time_extrapolation.h"
#include "time_integrator.hpp"
#include "../limex_tools.h"
//#include "../multi_thread_tools.h"
 
 
#ifdef UG_JSON
#include <nlohmann/json.hpp>
#endif
 
#undef LIMEX_MULTI_THREAD
 
namespace ug {
 
 
/*
template<class TI>
class TimeIntegratorThread : public TI
{
    TimeIntegratorThread(SmartPtr<grid_function_type> u1, number t1, ConstSmartPtr<grid_function_type> u0, number t0)
    {
 
    }
 
    // thre
    void apply()
    {
        TI::apply(u1, t1, u0, t0);
    }
 
 
};*/
 
static void MyPrintError(UGError &err)
{
    for(size_t i=0;i<err.num_msg();++i)
    {
        UG_LOG("MYERROR "<<i<<":"<<err.get_msg(i)<<std::endl);
        UG_LOG("     [at "<<err.get_file(i)<<
               ", line "<<err.get_line(i)<<"]\n");
    }
}
 
class ILimexRefiner
{
    virtual ~ILimexRefiner(){};
 
protected:
    SmartPtr<IRefiner> m_spRefiner;
};
class ILimexRefiner {…};


class ILimexCostStrategy
{
public:
    virtual ~ILimexCostStrategy(){}

    virtual void update_cost(std::vector<number> &costA, const std::vector<size_t> &vSteps, const size_t nstages) = 0;
};
class ILimexCostStrategy {…};
 

class LimexDefaultCost : public ILimexCostStrategy
{
public:
    LimexDefaultCost(){};
    void update_cost(std::vector<number> &m_costA, const std::vector<size_t> &m_vSteps, const size_t nstages)
    {
        UG_ASSERT(m_costA.size() >= nstages, "Huhh: Vectors match in size:" << m_costA.size() << "vs." << nstages);
 
        UG_LOG("A_0="<< m_vSteps[0] << std::endl);
        m_costA[0] = (1.0)*m_vSteps[0];
        for (size_t i=1; i<nstages; ++i)
        {
            m_costA[i] = m_costA[i-1] + (1.0)*m_vSteps[i];
            UG_LOG("A_i="<< m_vSteps[i] << std::endl);
        }
    }
    void update_cost(std::vector<number> &m_costA, const std::vector<size_t> &m_vSteps, const size_t nstages) {…}
};
class LimexDefaultCost : public ILimexCostStrategy {…};

class LimexNonlinearCost : public ILimexCostStrategy
{
public:
    LimexNonlinearCost() :
    m_cAssemble (1.0), m_cMatAdd(0.0), m_cSolution(1.0), m_useGamma(1)
    {}
 
    void update_cost(std::vector<number> &m_costA, const std::vector<size_t> &nSteps, const size_t nstages)
    {
        UG_ASSERT(m_costA.size() >= nstages, "Huhh: Vectors match in size:" << m_costA.size() << "vs." << nstages);
 
        // 3 assemblies (M0, J0, Gamma0)
        m_costA[0] = (2.0+m_useGamma)*m_cAssemble + nSteps[0]*((1.0+m_useGamma)*m_cMatAdd + m_cSolution);
        for (size_t i=1; i<=nstages; ++i)
        {
            // n-1 assemblies for Mk, 2n MatAdds, n solutions
            m_costA[i] = m_costA[i-1] + (nSteps[i]-1) * m_cAssemble + nSteps[i]*((1.0+m_useGamma)*m_cMatAdd + m_cSolution);
        }
    }
    void update_cost(std::vector<number> &m_costA, const std::vector<size_t> &nSteps, const size_t nstages) {…}
 
protected:
    double m_cAssemble;  
    double m_cMatAdd;   
    double m_cSolution;  
 
    int m_useGamma;
 
};
class LimexNonlinearCost : public ILimexCostStrategy {…};
 
 
class LimexTimeIntegratorConfig
 
{
public:
    LimexTimeIntegratorConfig()
    : m_nstages(2),
      m_tol(0.01),
      m_rhoSafety(0.8),
      m_sigmaReduction(0.5),
      m_greedyOrderIncrease(2.0),
      m_max_reductions(2),
      m_asymptotic_order(1000),
      m_useCachedMatrices(false),
      m_conservative(0)
      {}
 
    LimexTimeIntegratorConfig(unsigned int nstages)
    : m_nstages(nstages),
      m_tol(0.01),
      m_rhoSafety(0.8),
      m_sigmaReduction(0.5),
      m_greedyOrderIncrease(2.0),
      m_max_reductions(2),
      m_asymptotic_order(1000),
      m_useCachedMatrices(false),
      m_conservative(0)
        {}
 
 
protected:
    unsigned int m_nstages;                
    double m_tol;
    double m_rhoSafety;
    double m_sigmaReduction;
 
    double m_greedyOrderIncrease;
 
    size_t m_max_reductions;
    size_t m_asymptotic_order;                
 
    bool m_useCachedMatrices;
    unsigned int m_conservative; // stepping back?
 
#ifdef UG_JSON
        NLOHMANN_DEFINE_TYPE_INTRUSIVE(LimexTimeIntegratorConfig,
                        m_nstages, m_tol, m_rhoSafety,
                        m_sigmaReduction, m_greedyOrderIncrease, m_max_reductions, m_asymptotic_order,
                        m_useCachedMatrices, m_conservative)
#endif
public:
    std::string config_string() const {
    #ifdef UG_JSON
            try{
                nlohmann::json j(*this);
                    //  to_json(j, *this);
                return j.dump();
            } catch (...) {
                return std::string("EXCEPTION!!!");
            }
    #else
                    return std::string("LimexTimeIntegratorConfig::config_string()");
    #endif
        }
 
};
class LimexTimeIntegratorConfig {…};

template<class TDomain, class TAlgebra>
class LimexTimeIntegrator
: public INonlinearTimeIntegrator<TDomain, TAlgebra>,
  public DebugWritingObject<TAlgebra>,
  public LimexTimeIntegratorConfig // TODO: should become a 'has a'.
{
 
 
public:
        typedef TAlgebra algebra_type;
        typedef typename algebra_type::matrix_type matrix_type;
        typedef typename algebra_type::vector_type vector_type;
        typedef GridFunction<TDomain, TAlgebra> grid_function_type;
        typedef INonlinearTimeIntegrator<TDomain, TAlgebra> base_type;
        typedef typename base_type::solver_type solver_type;
 
        typedef IDomainDiscretization<algebra_type> domain_discretization_type;
        typedef LinearImplicitEuler<algebra_type> timestep_type;
        typedef AitkenNevilleTimex<vector_type> timex_type;
        typedef INonlinearTimeIntegrator<TDomain, TAlgebra> itime_integrator_type;
        typedef SimpleTimeIntegrator<TDomain, TAlgebra> time_integrator_type;
        typedef ISubDiagErrorEst<vector_type> error_estim_type;

        class ThreadData
        {
            //typedef boost::thread thread_type;
        public:
 
            ThreadData(SmartPtr<timestep_type> spTimeStep)
            : m_stepper(spTimeStep)
            {}
 
            ThreadData(){}
            SmartPtr<timestep_type> get_time_stepper()
            { return m_stepper; }
 
 
            void set_solver(SmartPtr<solver_type> solver)
            { m_solver = solver;}
 
            SmartPtr<solver_type> get_solver()
            { return m_solver;}
 
            void set_error(int e)
            { m_error=e; }
 
 
            void get_error()
            { /*return m_error;*/ }
 
 
            void set_solution(SmartPtr<grid_function_type> sol)
            { m_sol = sol;}
 
            SmartPtr<grid_function_type> get_solution()
            { return m_sol;}
 
            void set_derivative(SmartPtr<grid_function_type> sol)
            { m_dot = sol;}
 
            SmartPtr<grid_function_type> get_derivative()
            { return m_dot;}
 
        protected:
             // includes time step series
            SmartPtr<timestep_type> m_stepper;
            SmartPtr<solver_type> m_solver;
 
            SmartPtr<grid_function_type> m_sol;
            SmartPtr<grid_function_type> m_dot;
            int m_error;
        };
        class ThreadData {…};
 
        typedef std::vector<SmartPtr<ThreadData> > thread_vector_type;
 
public:

    void set_debug_for_timestepper(SmartPtr<IDebugWriter<algebra_type> > spDebugWriter)
    {
        for(size_t i=0; i<m_vThreadData.size(); ++i)
        {
                m_vThreadData[i].get_time_stepper()->set_debug(spDebugWriter);
        }
        UG_LOG("set_debug:" << m_vThreadData.size());
    }
    void set_debug_for_timestepper(SmartPtr<IDebugWriter<algebra_type> > spDebugWriter) {…}
 
    //using VectorDebugWritingObject<vector_type>::set_debug;
protected:
    //using VectorDebugWritingObject<vector_type>::write_debug;
 
 
 
public:
        LimexTimeIntegrator(int nstages):
          LimexTimeIntegratorConfig(nstages),
          m_gamma(m_nstages+1),
          m_costA(m_nstages+1),
          m_monitor(((m_nstages)*(m_nstages))), // TODO: wasting memory here!
          m_workload(m_nstages),
          m_lambda(m_nstages), 
          m_num_reductions(m_nstages, 0),
          m_consistency_error(m_nstages),
          m_spCostStrategy(make_sp<LimexDefaultCost>(new LimexDefaultCost())),
          m_spBanachSpace(new AlgebraicSpace<grid_function_type>()),              // default algebraic space
          m_bInterrupt(false),
          m_limex_step(1)
        {
            m_vThreadData.reserve(m_nstages);
            m_vSteps.reserve(m_nstages);
            init_gamma(); // init exponents (i.e. k+1, k, 2k+1, ...)
        }
 
 

        void set_tolerance(double tol) { m_tol = tol;}
        void set_stepsize_safety_factor(double rho) { m_rhoSafety = rho;}
        void set_stepsize_reduction_factor(double sigma) { m_sigmaReduction = sigma;}
        void set_stepsize_greedy_order_factor(double sigma) { m_greedyOrderIncrease = sigma;}
 
        void set_max_reductions(size_t nred) { m_max_reductions = nred;}
        void set_asymptotic_order(size_t q) { m_asymptotic_order = q;}
 
        void set_start_step(size_t step){m_limex_step=step;}

        void add_error_estimator(SmartPtr<error_estim_type> spErrorEstim)
        { m_spErrorEstimator = spErrorEstim; }
        void add_error_estimator(SmartPtr<error_estim_type> spErrorEstim) {…}

        void add_stage_base(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD, SmartPtr<domain_discretization_type> spGamma=SPNULL)
        {
            UG_ASSERT(m_vThreadData.size() == m_vSteps.size(), "ERROR: m_vThreadData and m_vSteps differ in size!");
 
            UG_ASSERT(m_vThreadData.empty() || m_vSteps.back()<nsteps, "ERROR: Sequence of steps must be increasing." );
 
            // all entries have been set
            if (m_vThreadData.size() == m_nstages)
            { return; }
 
 
            // a) set number of steps
            m_vSteps.push_back(nsteps);
 
            // b) set time-stepper
            SmartPtr<timestep_type> limexStepSingleton;
#ifndef LIMEX_MULTI_THREAD
 
            if(m_vThreadData.size()>0) {
                // re-use time-stepper (if applicable)
                limexStepSingleton = m_vThreadData.back().get_time_stepper();
            }
            else
            {
                // create time-stepper
#endif
                // for mult-threading, each info object has own time-stepper
                if (spGamma.invalid()) {
                    limexStepSingleton = make_sp(new timestep_type(spDD));
                } else {
                    limexStepSingleton = make_sp(new timestep_type(spDD, spDD, spGamma));
                }
                UG_ASSERT(limexStepSingleton.valid(), "Huhh: Invalid pointer")
#ifndef LIMEX_MULTI_THREAD
            }
#endif
            // propagate debug info
            limexStepSingleton->set_debug(VectorDebugWritingObject<vector_type>::vector_debug_writer());
            m_vThreadData.push_back(ThreadData(limexStepSingleton));
 
            // c) set solver
            UG_ASSERT(solver.valid(), "Huhh: Need to supply solver!");
            m_vThreadData.back().set_solver(solver);
            UG_ASSERT(m_vThreadData.back().get_solver().valid(), "Huhh: Need to supply solver!");
        }
        void add_stage_base(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD, SmartPtr<domain_discretization_type> spGamma=SPNULL) {…}
 
        void add_stage(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD)
        { add_stage_base(nsteps, solver, spDD); }
 
        void add_stage_ext(size_t nsteps, SmartPtr<solver_type> solver, SmartPtr<domain_discretization_type> spDD, SmartPtr<domain_discretization_type> spGamma)
        { add_stage_base(nsteps, solver, spDD, spGamma); }

        void add_stage(size_t i, size_t nsteps, SmartPtr<domain_discretization_type> spDD, SmartPtr<solver_type> solver)
        {
            UG_LOG("WARNING: add_stage(i, nsteps ,...) is deprecated. Please use 'add_stage(nsteps ,...) instead!'");
            add_stage(nsteps, solver, spDD);
        }
        void add_stage(size_t i, size_t nsteps, SmartPtr<domain_discretization_type> spDD, SmartPtr<solver_type> solver) {…}

        SmartPtr<timestep_type> get_time_stepper(size_t i)
        {
            UG_ASSERT(i<m_vThreadData.size(), "Huhh: Invalid entry");
            return m_vThreadData[i].get_time_stepper();
        }
 
 
 
 
protected:
        void init_integrator_threads(ConstSmartPtr<grid_function_type> u);

        int apply_integrator_threads(number dtcurr, ConstSmartPtr<grid_function_type> u0, number t0, size_t nstages);

        void join_integrator_threads();

        void update_integrator_threads(ConstSmartPtr<grid_function_type> ucommon, number t);

        void dispose_integrator_threads();
 
public:
        bool apply(SmartPtr<grid_function_type> u, number t1, ConstSmartPtr<grid_function_type> u0, number t0);
 
        number  get_cost(size_t i) { return m_costA[i]; }
        number  get_gamma(size_t i) { return m_gamma[i]; }
        number  get_workload(size_t i) { return m_workload[i]; }
 
 
protected:
        number& monitor(size_t k, size_t q) {
          UG_ASSERT(k<m_nstages, "Huhh: k mismatch");
          UG_ASSERT(q<m_nstages, "Huhh: q mismatch");
          return m_monitor[k+(m_nstages)*q]; 
        }

        void init_gamma()
        {
            for (size_t k=0; k<=m_nstages; ++k)
            {
                m_gamma[k] = 2.0+k;
            }
        }
        void init_gamma() {…}

        //  (depends on m_vSteps, which must have been initialized!)
        void update_cost()
        {
            m_spCostStrategy->update_cost(m_costA, m_vSteps, m_nstages);
        }
        void update_cost() {…}

        // (depends on cost, which must have been initialized!)
        void update_monitor()
        {
          UG_ASSERT(m_costA.size()>= m_nstages, "Cost vector too small!")
          for (size_t k=0; k<m_nstages-1; ++k)
          {
            UG_LOG("k= "<< k<< ", A[k]=" << m_costA[k] << ", gamma[k]=" << m_gamma[k] << "\t");
            for (size_t q=0; q<m_nstages-1; ++q)
              {
            // Deuflhard: Order and stepsize, ... eq. (3.7)
            double gamma = (m_costA[k+1] - m_costA[0] + 1.0)/(m_costA[q+1] - m_costA[0] + 1.0);
            double alpha = pow(m_tol, gamma);
            
            // for fixed order q, the monitor indicates the performance penalty compared to a strategy using k stages only
            // cf. eq. (4.6)
            monitor(k,q) = pow(alpha/(m_tol * m_rhoSafety), 1.0/m_gamma[k]);
            UG_LOG(monitor(k,q) << "[" << pow(alpha/(m_tol), 1.0/m_gamma[k]) <<"," << gamma<< ","<< m_costA[k+1] << "," << m_costA[q+1]<< ","<< alpha << "]" << "\t");
            // UG_LOG(  << "\t");
            
              }
            UG_LOG(std::endl);
          }
 
        }
        void update_monitor() {…}
 
        // Find row k=1, ..., ntest-1 minimizing (estimated) error eps[kmin]
        // Also: predict column q with minimal workload W_{k+1,k} = A_{k+1} * lambda_{k+1}
        size_t find_optimal_solution(const std::vector<number>& eps, size_t ntest, /*size_t &kf,*/ size_t &qpred);
 
public:
        void set_time_derivative(SmartPtr<grid_function_type> udot) {m_spDtSol = udot;}

        SmartPtr<grid_function_type> get_time_derivative() {return m_spDtSol;}

        bool has_time_derivative() {return m_spDtSol!=SPNULL;}
 
 
        void enable_matrix_cache() { m_useCachedMatrices = true; }   
        void disable_matrix_cache() { m_useCachedMatrices = false; }    
 
        void select_cost_strategy(SmartPtr<ILimexCostStrategy> cost)
        { m_spCostStrategy = cost;}
 

        void set_space(SmartPtr<IGridFunctionSpace<grid_function_type> > spSpace)
        { 
          m_spBanachSpace = spSpace;
          UG_LOG("set_space:" << m_spBanachSpace->config_string()); 
        }
        void set_space(SmartPtr<IGridFunctionSpace<grid_function_type> > spSpace) {…}

        void interrupt() {m_bInterrupt = true;}
 
        void set_conservative(bool c)
        { m_conservative = (c) ? 1 : 0; }
 
        std::string config_string() const
        { return LimexTimeIntegratorConfig::config_string(); }
 
protected:
 
        SmartPtr<error_estim_type> m_spErrorEstimator;     // (smart ptr for) error estimator
 
        std::vector<size_t> m_vSteps;           
        std::vector<ThreadData> m_vThreadData; 
 
        std::vector<number> m_gamma;         
 
        std::vector<number> m_costA;         
        std::vector<number> m_monitor;         
 
        std::vector<number> m_workload;
        std::vector<number> m_lambda;
 
        std::vector<size_t> m_num_reductions;       
 
        std::vector<number> m_consistency_error; 
 
        SmartPtr<grid_function_type> m_spDtSol;    
 
        SmartPtr<ILimexCostStrategy> m_spCostStrategy;

        SmartPtr<IGridFunctionSpace<grid_function_type> > m_spBanachSpace;
 
        bool m_bInterrupt;
        int m_limex_step;                       
 
 
};
 

template<class TDomain, class TAlgebra>
void LimexTimeIntegrator<TDomain,TAlgebra>::init_integrator_threads(ConstSmartPtr<grid_function_type> u)
{
    PROFILE_FUNC_GROUP("limex");
    const int nstages = m_vThreadData.size()-1;
    for (int i=nstages; i>=0; --i)
    {
        m_vThreadData[i].set_solution(u->clone());
        m_vThreadData[i].set_derivative(u->clone());
        m_vThreadData[i].get_time_stepper()->set_matrix_cache(m_useCachedMatrices);
    }
}
void LimexTimeIntegrator<TDomain,TAlgebra>::init_integrator_threads(ConstSmartPtr<grid_function_type> u) {…}

template<class TDomain, class TAlgebra>
void LimexTimeIntegrator<TDomain,TAlgebra>::dispose_integrator_threads()
{
    PROFILE_FUNC_GROUP("limex");
    const int nstages = m_vThreadData.size()-1;
    for (int i=nstages; i>=0; --i)
    {
        m_vThreadData[i].set_solution(SPNULL);
        m_vThreadData[i].set_derivative(SPNULL);
        // m_vThreadData[i].get_time_stepper()->set_matrix_cache(m_useCachedMatrices);
    }
}
void LimexTimeIntegrator<TDomain,TAlgebra>::dispose_integrator_threads() {…}
 
 
// create (& execute) threads
/*boost::thread_group g;
        typename thread_vector_type::reverse_iterator rit=m_vThreadData.rbegin();
        for (rit++; rit!= m_vThreadData.rend(); ++rit)
        {
 
            boost::thread *t =new boost::thread(boost::bind(&ThreadSafeTimeIntegrator::apply, *rit));
            //g.add_thread(t);
 
            g.create_thread(boost::bind(&ThreadSafeTimeIntegrator::apply, *rit));
 
        }*/
 
 
/* TODO: PARALLEL execution?*/
template<class TDomain, class TAlgebra>
int LimexTimeIntegrator<TDomain,TAlgebra>::apply_integrator_threads(number dtcurr, ConstSmartPtr<grid_function_type> u0, number t0, size_t nstages)
{
    PROFILE_FUNC_GROUP("limex");
    update_cost();       // compute cost A_i (alternative: measure times?)
    update_monitor(); // convergence monitor
 
    /*
            int tn = omp_get_thread_num();
            int nt = omp_get_num_threads();
            omp_set_num_threads(nstages);
     */
    int error = 0;
    //const int nstages = m_vThreadData.size()-1;
    //  #pragma omp for private(i) // shared (nstages, u1) schedule(static)
    for (int i=0; i<=(int)nstages; ++i)
    {
        /*
                std::cerr << "I am " << tn << " of " << nt << " ("<< i<< "/" << nstages<<")!" << std::endl;
                UGMultiThreadEnvironment mt_env;
         */
 
        // copy data to private structure (one-to-many)
        //m_vThreadData[i].set_solution(u1->clone());
 
        // switch to "child" comm
        // mt_env.begin();
 
        // integrate (t0, t0+dtcurr)
        time_integrator_type integrator(m_vThreadData[i].get_time_stepper());
        integrator.set_time_step(dtcurr/m_vSteps[i]);
        integrator.set_dt_min(dtcurr/m_vSteps[i]);
        integrator.set_dt_max(dtcurr/m_vSteps[i]);
        integrator.set_reduction_factor(0.0);                 // quit immediately, if step fails
        integrator.set_solver(m_vThreadData[i].get_solver());
        integrator.set_derivative(m_vThreadData[i].get_derivative());
        
        if(this->debug_writer_valid())
        {
            integrator.set_debug(this->debug_writer());
            char debug_name_ext[16]; snprintf(debug_name_ext, 16, "%04d", i);
            this->enter_debug_writer_section(std::string("Stage_") + debug_name_ext);
        }
 
        UG_ASSERT(m_spBanachSpace.valid(), "Huhh: Need valid (default) banach space");
        integrator.set_banach_space(m_spBanachSpace);
        UG_LOG("Set space:" << m_spBanachSpace->config_string());
 
        bool exec = true;
        try
        {
            exec = integrator.apply(m_vThreadData[i].get_solution(), t0+dtcurr, u0, t0);
            m_consistency_error[i] = integrator.get_consistency_error();
        }
        catch(ug::UGError& err)
        {
 
            exec = false;
            error += (1 << i);
            UG_LOGN("Step "<< i<< " failed on stage " << i << ": " << err.get_msg());
            MyPrintError(err);
 
        }
        
        this->leave_debug_writer_section();
 
        if (!exec)
        {
            // Additional actions at failure
            UG_LOGN("Step "<< i<< " failed on stage " << i << ": reason not clear!" );
        }
 
        // switch to "parent" comm
        //mt_env.end();
    } /*for all stages loop*/
 
 
 
    return error;
}
int LimexTimeIntegrator<TDomain,TAlgebra>::apply_integrator_threads(number dtcurr, ConstSmartPtr<grid_function_type> u0, number t0, size_t nstages) {…}
 
 
template<class TDomain, class TAlgebra>
void LimexTimeIntegrator<TDomain,TAlgebra>::join_integrator_threads()
{
    // join all threads
    // g.join_all();
    const int nstages = m_vThreadData.size()-1;
    for (int i=nstages; i>=0; --i)
    {
        m_vThreadData[i].get_time_stepper()->invalidate();
    }
}
void LimexTimeIntegrator<TDomain,TAlgebra>::join_integrator_threads() {…}
 
 
template<class TDomain, class TAlgebra>
void LimexTimeIntegrator<TDomain,TAlgebra>::update_integrator_threads(ConstSmartPtr<grid_function_type> ucommon, number t)
{
    const int nstages = m_vThreadData.size()-1;
    for (int i=nstages; i>=0; --i)
    {
        UG_ASSERT(m_vThreadData[i].get_solution()->size()==ucommon->size(), "LIMEX: Vectors must match in size!")
        *m_vThreadData[i].get_solution() = *ucommon;
    }
}
void LimexTimeIntegrator<TDomain,TAlgebra>::update_integrator_threads(ConstSmartPtr<grid_function_type> ucommon, number t) {…}
 
template<class TDomain, class TAlgebra>
size_t LimexTimeIntegrator<TDomain,TAlgebra>::
find_optimal_solution(const std::vector<number>& eps, size_t ntest, /*size_t &kf,*/ size_t &qpred)
{
 
    const size_t qold=qpred;
 
    size_t jbest = 1;
    qpred = 1;
 
    size_t j=1;
    size_t k=j-1;
 
    m_lambda[k] = pow(m_rhoSafety*m_tol/eps[j], 1.0/m_gamma[k]);   // 1/epsilon(k)
    m_workload[k] = m_costA[j]/m_lambda[k];
    UG_LOG("j=" << j << ": eps=" << eps[j]  << ", lambda(j)=" <<m_lambda[k]  << ", epsilon(j)=" <<1.0/m_lambda[k] << "<= alpha(k, qcurr)=" << monitor(k,qold-1) << "< alpha(k, qcurr+1)=" << monitor(k,qold) <<", A="<< m_costA[j] << ", W="<< m_workload[k] <<std::endl);
 
    for (j=2; j<ntest; ++j)
    {
      k = j-1;
      m_lambda[k] = pow(m_rhoSafety*m_tol/eps[j], 1.0/m_gamma[k]);
      m_workload[k] = m_costA[j]/m_lambda[k];
      UG_LOG("j=" << j << ": eps=" << eps[j]  << ", lambda(j)=" <<m_lambda[k]  << ", epsilon(j)=" <<1.0/m_lambda[k] << "<= alpha(k, qcurr)=" << monitor(k,qold-1) << "< alpha(k, qcurr+1)=" << monitor(k,qold) <<", A="<< m_costA[j] << ", W="<< m_workload[k] <<std::endl);
 
 
      // TODO: Convergence monitor
 
      qpred = (m_workload[qpred-1] > m_workload[k]) ? j : qpred;
      jbest = (eps[jbest] > eps [j]) ? j : jbest;
    }
 
    return jbest;
}
 
template<class TDomain, class TAlgebra>
bool LimexTimeIntegrator<TDomain,TAlgebra>::
apply(SmartPtr<grid_function_type> u, number t1, ConstSmartPtr<grid_function_type> u0, number t0)
{
    PROFILE_FUNC_GROUP("limex");
#ifdef UG_OPENMP
    // create multi-threading environment
    //int nt = std::min(omp_get_max_threads(), m_nstages);
 
#endif
 
    // NOTE: we use u as common storage for future (and intermediate) solution(s)
    if (u.get() != u0.get())    // only copy if not already identical, otherwise: PST_UNDEFINED!
        *u = *u0;
 
    // initialize integrator threads
    // (w/ solutions)
    init_integrator_threads(u0);
 
 
    // write_debug
    for (unsigned int i=0; i<m_vThreadData.size(); ++i)
    {
        std::ostringstream ossName;
        ossName << std::setfill('0') << std::setw(4);
        ossName << "Limex_Init_iter" << 0 << "_stage" << i;
        this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
    }
 
    number t = t0;
    double dtcurr = ITimeIntegrator<TDomain, TAlgebra>::get_time_step();
 
    size_t kmax = m_vThreadData.size();        // maximum number of stages
    size_t qpred = kmax-1;                              // predicted optimal order
    size_t qcurr = qpred;                               // current order
 
    // for Gustafsson/lundh/Soederlind type PID controller
    /*size_t qlast = 0;
    double dtlast = 0.0;
    std::vector<double> epslast(kmax, m_rhoSafety*m_tol);
*/
 
    // time integration loop
    SmartPtr<grid_function_type> ubest = SPNULL;
    size_t limex_total = 1;
    size_t limex_success = 0;
    size_t ntest;    
    size_t jbest;
 
    m_bInterrupt = false;
    //bool bProbation = false;
    bool bAsymptoticReduction = false;
 
    const size_t nSwitchHistory=16;
    const size_t nSwitchLookBack=5;
    int vSwitchHistory[nSwitchHistory] ={ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
 
    timex_type timex(m_vSteps);
    while ((t < t1) && ((t1-t) > base_type::m_precisionBound))
    {
        int err = 0;
 
        //UG_DLOG(LIB_LIMEX, 5, "+++ LimexTimestep +++" << limex_step << "\n");
        UG_LOG("+++ LimexTimestep +++" << m_limex_step << "\n");
 
 
        // save time stamp for limex step start
        Stopwatch stopwatch;
        stopwatch.start();
 
 
        // determine step size
        number dt = std::min(dtcurr, t1-t);
        UG_COND_THROW(dt < base_type::get_dt_min(), "Time step size below minimum. ABORTING! dt=" << dt << "; dt_min=" << base_type::get_dt_min() << "\n");
 
 
        // Notify init observers. (NOTE: u = u(t))¸
        itime_integrator_type::notify_init_step(u, m_limex_step, t, dt);
 
 
        // determine number of stages to investigate
        qcurr = qpred;
        ntest = std::min(kmax, qcurr+1);
        UG_LOG("ntest="<< ntest << std::endl);
 
        // checks
        UG_ASSERT(m_vSteps.size() >= ntest, "Huhh: sizes do not match: " << m_vSteps.size() << "<"<<ntest);
        UG_ASSERT(m_vThreadData.size() >= ntest, "Huhh: sizes do not match: " << m_vThreadData.size() << "< "<< ntest);

        // PARALLEL EXECUTION: BEGIN
 
        // write_debug
        for (size_t i=0; i<ntest; ++i)
        {
            std::ostringstream ossName;
            ossName << std::setfill('0') << std::setw(4);
            ossName << "Limex_BeforeSerial_iter" << m_limex_step << "_stage" << i << "_total" << limex_total;
            this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
        }
        if(this->debug_writer_valid())
        {
            char debug_step_id_ext[16]; snprintf(debug_step_id_ext, 16, "%04d", m_limex_step);
            char debug_total_id_ext[16]; snprintf(debug_total_id_ext, 16, "%04d", (int) limex_total);
            this->enter_debug_writer_section(std::string("LimexTimeIntegrator_iter") + debug_step_id_ext + "_total" + debug_total_id_ext);
        }
 
        // integrate: t -> t+dt
        err = apply_integrator_threads(dt, u, t, ntest-1);
        
        this->leave_debug_writer_section();
 
        // write_debug
        for (size_t i=0; i<ntest; ++i)
        {
            std::ostringstream ossName;
            ossName << std::setfill('0') << std::setw(4);
            ossName << "Limex_AfterSerial_iter" << m_limex_step << "_stage" << i << "_total" << limex_total;
            this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
        }
 
        join_integrator_threads();
        // PARALLEL EXECUTION: END
 

        // SERIAL EXECUTION: BEGIN
 
        // sanity checks
        UG_ASSERT(m_spErrorEstimator.valid(), "Huhh: Invalid Error estimator?");
 
        double epsmin = 0.0;
 
        bool limexConverged = false;
        if (err==0)
        {   // Compute extrapolation at t+dtcurr (SERIAL)
            // Consistency check.
            for (unsigned int i=1; i<ntest; ++i)
            {
                UG_LOG("Checking consistency:" <<m_consistency_error[i] <<"/" << m_consistency_error[i-1] << "="<< m_consistency_error[i]/m_consistency_error[i-1] << std::endl);
            }
 
            // Extrapolation.
            timex.set_error_estimate(m_spErrorEstimator);
            for (unsigned int i=0; i<ntest; ++i)
            {
                UG_ASSERT(m_vThreadData[i].get_solution().valid(), "Huhh: no valid solution?");
                timex.set_solution(m_vThreadData[i].get_solution(), i);
            }
            timex.apply(ntest);
 
            // write_debug
            for (size_t i=0; i<ntest; ++i)
            {
                std::ostringstream ossName;
                ossName << std::setfill('0') << std::setw(4);
                ossName << "Limex_Extrapolates_iter" << m_limex_step << "_stage" << i << "_total" << limex_total;
                this->write_debug(*m_vThreadData[i].get_solution(), ossName.str().c_str());
            }
            limex_total++;
 
            // Obtain sub-diagonal error estimates.
            const std::vector<number>& eps = timex.get_error_estimates();
            UG_ASSERT(ntest<=eps.size(), "Huhh: Not enough solutions?");
 
            // select optimal solution (w.r.t error) AND
            // predict optimal order (w.r.t. workload) for next step
            jbest = find_optimal_solution(eps, ntest, qpred);
            UG_ASSERT(jbest < ntest, "Huhh: Not enough solutions?");
 
            // best solution
            ubest  = timex.get_solution(jbest-m_conservative).template cast_dynamic<grid_function_type>();
            epsmin = eps[jbest];
 
            // check for convergence
            limexConverged = (epsmin <= m_tol) ;
 
 
            if (limex_success>3)
            {
 
                vSwitchHistory[m_limex_step%nSwitchHistory] = (qpred - qcurr);
                UG_DLOG(LIB_LIMEX, 5, "LIMEX-ASYMPTOTIC-ORDER switch:  = " <<  (qpred - qcurr)<< std::endl);
 
                size_t nSwitches=0;
                for (int s=nSwitchLookBack-1; s>=0; s--)
                {
                    nSwitches += std::abs(vSwitchHistory[(m_limex_step-s)%nSwitchHistory]);
                    UG_DLOG(LIB_LIMEX, 6, "LIMEX-ASYMPTOTIC-ORDER: s[" << s<< "] = " <<  vSwitchHistory[(m_limex_step-s)%nSwitchHistory] << std::endl);
                }
                UG_DLOG(LIB_LIMEX, 5,"LIMEX-ASYMPTOTIC-ORDER: nSwitches = " <<  nSwitches << std::endl);
 
 
            /*
            if (bProbation)
            {
 
                if ((qpred < qcurr) || (!limexConverged))   // consecutive (follow-up) order decrease
                {
                    bProbation = true;
                    m_num_reductions[0]++;
                    UG_LOG("LIMEX-ASYMPTOTIC-ORDER: Decrease on parole detected: "<< qcurr << " -> " << qpred << "("<< m_num_reductions[0]<<")"<<std::endl);
 
                }
                else
                {
                    // reset all counters
                    bProbation = false;
                    UG_LOG("LIMEX-ASYMPTOTIC-ORDER: Probation off!"<<std::endl);
 
                }
            }
            else
            {
                if ((qpred < qcurr) || (!limexConverged))// 1st order decrease
                {
                    bProbation = true;
                    m_num_reductions[0]++;
                    UG_LOG("LIMEX-ASYMPTOTIC-ORDER:  Decrease detected: "<< qcurr << " -> " << qpred << "("<< m_num_reductions[0]<<")"<<std::endl);
                }
                else
                {
                    m_num_reductions[0] = 0;
                    UG_LOG("LIMEX-ASYMPTOTIC-ORDER: I am totally free!"<<std::endl);
                }
 
 
            }
*/
 
            // bAsymptoticReduction = (m_num_reductions[0] >= m_max_reductions) || bAsymptoticReduction;
 
            // bAsymptoticReduction = (nSwitches >= m_max_reductions) || bAsymptoticReduction;
 
            if (nSwitches >= m_max_reductions) {
 
 
                m_asymptotic_order = std::min<size_t>(m_asymptotic_order-1, 2);
                bAsymptoticReduction = true;
 
                for (int s=nSwitchLookBack-1; s>=0; s--)
                { vSwitchHistory[(m_limex_step-s)%nSwitchHistory]=0; }
 
 
            }
 
 
 
            // asymptotic order reduction
            //if (m_num_reductions[0] >= m_max_reductions)
            if  (bAsymptoticReduction)
            {
                kmax = (kmax >= m_asymptotic_order ) ? m_asymptotic_order : kmax;
                qpred = kmax - 1;
                UG_DLOG(LIB_LIMEX, 5, "LIMEX-ASYMPTOTIC-ORDER: Reduction: "<< qpred );
                UG_DLOG(LIB_LIMEX, 5, "(kmax=" << kmax << ", asymptotic"<<  m_asymptotic_order <<") after "<<m_max_reductions<< std::endl);
            }
 
            } // (limex_success>3)
 
            /*
            // adjust for order bound history
            if ((m_num_reductions[qpred] > m_max_reductions) && (qpred > 1))
            {
                    UG_LOG("prohibited  q:"<< qpred << "(" << m_num_reductions[qpred] << ")"<<std::endl);
                    //qpred--;
                    qpred = kmax = 2;
            } else
            {
                UG_LOG("keeping  q:"<< qpred << "(" << m_num_reductions[qpred] << ")"<<std::endl);
            }
            */
 
            //double pid = 1.0;
 
            // constant order?
            /*if (qcurr == qlast)
            {
                double dtRatio =dtcurr/dtlast;
                // Dd, Band 3 (DE), p. 360
                int k = ntest-1;
                while (k--) {
                    double epsRatio =eps[k+1]/epslast[k+1];
 
                    double qopt = log(epsRatio) / log(dtRatio);  // take dtcurr here, as dt may be smaller
 
                    UG_LOG("effective p["<< k << "]="<< qopt-1.0 <<","<< eps[k+1] <<","<< epslast[k+1] <<","  <<  epsRatio << "("<<  log(epsRatio)  << ","<<  pow(epsRatio, 1.0/m_gamma[k])  <<")"<< dtcurr <<","<< dtlast   << "," << dtRatio << "," <<log(dtRatio) << std::endl);
                }
 
 
                double epsRatio = eps[qcurr]/epslast[qcurr];
                pid = dtRatio/pow(epsRatio, 1.0/m_gamma[qcurr-1]);
 
                UG_LOG("pid=" << pid << std::endl);
 
            }*/
 
 
            // select (predicted) order for next step//double dtpred = dtcurr*std::min(m_lambda[qpred-1], itime_integrator_type::get_increase_factor());
            double dtpred = dtcurr*std::min(m_lambda[qpred-1], itime_integrator_type::get_increase_factor());
            //double dtpred = dtcurr*m_lambda[qpred-1];
            UG_LOG("+++++\nget_increase_factor() gives "<<itime_integrator_type::get_increase_factor()<<" \n+++++++")
            UG_LOG("koptim=\t" << jbest << ",\t eps(k)=" << epsmin << ",\t q=\t" << qpred<< "("<<  ntest << "), lambda(q)=" << m_lambda[qpred-1] << ", alpha(q-1,q)=" << monitor(qpred-1, qpred) << "dt(q)=" << dtpred<< std::endl);
 
            // EXTENSIONS: convergence model
            if (limexConverged)
            {
                limex_success++;
 
                // a) aim for order increase in next step
              if ((qpred+1==ntest)  /* increase by one possible? */
                  //    && (m_lambda[qpred-1]>1.0)
                  // && (m_num_reductions[qpred+1] <= m_max_reductions)
                  && (kmax>ntest)) /* still below max? */
                {
                    const double alpha = monitor(qpred-1, qpred);
                    UG_LOG("CHECKING for order increase: "<< m_costA[qpred] << "*" << alpha << ">" << m_costA[qpred+1]);
                    // check, whether further increase could still be efficient
                    if (m_costA[qpred] * alpha > m_costA[qpred+1])
                    {
                        qpred++;                // go for higher order
                        if (m_greedyOrderIncrease >0.0) {
                          dtpred *= (1.0-m_greedyOrderIncrease) + m_greedyOrderIncrease*alpha;      // & adapt time step  // TODO: check required!
                        }
                        UG_LOG("... yes.\n")
 
                        // update history
                        vSwitchHistory[m_limex_step%nSwitchHistory] = (qpred - qcurr);
                        UG_LOG("LIMEX-ASYMPTOTIC-ORDER switch update:  = " <<  (qpred - qcurr)<< std::endl);
 
                    } else {
                        UG_LOG("... nope.\n")
                    }
 
                }
 
 
                // b) monitor convergence (a-priori check!)
 
            }
 
            // parameters for subsequent step
            dtcurr = std::min(dtpred, itime_integrator_type::get_dt_max());
 
 
        }
        else
        {
            // solver failed -> cut time step //
            number base_dtmin=base_type::get_dt_min();
            if(dtcurr <= base_dtmin)
                    base_type::set_dt_min(base_dtmin*m_sigmaReduction);
            
            dtcurr=std::max(base_type::get_dt_min(), dtcurr*m_sigmaReduction);
            
            //dtcurr=std:max(base_type::get_dt_min(), dtcurr);
        }
 
        // output compute time for Limex step
        number watchTime = stopwatch.ms()/1000.0;
        UG_LOGN("Time: " << std::setprecision(3) << watchTime << "s");
 
        if ((err==0) && limexConverged)
        {
            // ACCEPT time step
            UG_LOG("+++ LimexTimestep +++" << m_limex_step << " ACCEPTED"<< std::endl);
            UG_LOG("               :\t time \t dt (success) \t dt (pred) \tq=\t order (curr)" << qcurr+1 << std::endl);
            UG_LOG("LIMEX-ACCEPTING:\t" << t <<"\t"<< dt << "\t" << dtcurr << "\tq=\t" << qcurr+1 << std::endl);
 
            // update PID controller
            /*qlast = qcurr;
            epslast =  timex.get_error_estimates();  // double check ???
            dtlast = dt;*/
            // Compute time derivatives (by extrapolation)
            if (this->has_time_derivative())
            {
                UG_LOG("Computing derivative" << std::endl);
                grid_function_type &udot = *get_time_derivative();
 
                for (size_t i=0; i<=jbest; ++i)
                {
                    timex.set_solution(m_vThreadData[i].get_derivative(), i);
                }
                timex.apply(jbest+1, false); // do not compute error
 
                udot = *timex.get_solution(jbest).template cast_dynamic<grid_function_type>();
 
                std::ostringstream ossName;
                ossName << std::setfill('0');
                ossName << "Limex_Derivative_iter" << m_limex_step << "_total" << limex_total;
                this->write_debug(udot, ossName.str().c_str());
            }
 
 
            // Copy best solution.
            UG_ASSERT(ubest.valid(), "Huhh: Invalid error estimate?");
            *u = *ubest;
            t += dt;
 
            // Initiate post process.
            itime_integrator_type::notify_finalize_step(u, m_limex_step++, t, dt);
 
            // make sure that all threads continue
            // with identical initial value u(t)
            // update_integrator_threads(ubest, t);
 
 
            // working on last row => increase order
            //if (ntest == q+1) ntest++;
        }
        else
        {
            // DISCARD time step
            UG_LOG("+++ LimexTimestep +++" << m_limex_step << " FAILED" << std::endl);
            UG_LOG("               :\t time \t dt (failed) \t dt (curr) \teps=\t" << epsmin <<"\t(tol="<<m_tol<<")\terr="<<err<< std::endl);
            UG_LOG("LIMEX-REJECTING:\t" << t <<"\t"<< dt << "\t" << dtcurr <<std::endl);
 
            itime_integrator_type::notify_rewind_step(ubest, m_limex_step, t+dt, dt);
 
        }
 
        // interrupt if interruption flag set
        if (m_bInterrupt)
        {
            UG_LOGN("Limex interrupted by external command.");
            break;
        }
 
        // SERIAL EXECUTION: END
        update_integrator_threads(u, t);
 
 
        // SOLVE
 
        // ESTIMATE
 
        // MARK
 
        // REFINE
 
 
    } // time integration loop
 
    dispose_integrator_threads();
 
    return true;
} // apply
bool LimexTimeIntegrator<TDomain,TAlgebra>:: {…}
 
 
} // namespace ug
 
#endif /* TIME_INTEGRATOR_HPP_ */
        void disable_matrix_cache() { m_useCachedMatrices = false; }     {…}
class LimexTimeIntegrator {…};