time_integrator_impl.hpp

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: Copyright (c) 2014:  G-CSC, Goethe University Frankfurt
 * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
 *
 * Author: Arne Naegel, Andreas Kreienbuehl
 *
 * 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 __H__LIMEX__TIME_INTEGRATOR_IMPL_HPP__
#define __H__LIMEX__TIME_INTEGRATOR_IMPL_HPP__
 
namespace ug {
 
template<typename TDomain, typename TAlgebra>
bool LinearTimeIntegrator<TDomain, TAlgebra>::apply
(
  SmartPtr<grid_function_type> u1,
  number t1,
  ConstSmartPtr<grid_function_type> u0,
  number t0
)
{
 
  LIMEX_PROFILE_FUNC()
  UG_COND_THROW(!base_type::m_spLinearSolver.valid(), "Linear solver invalid");
 
  // short-cuts
  GridLevel const &gl = u0->grid_level();
  time_disc_type &tdisc = *tdisc_dep_type::m_spTimeDisc;
 
  // create solution vector & right hand side
  SmartPtr<typename base_type::vector_type> uold= u0->clone();
  SmartPtr<typename base_type::vector_type> b= u0->clone_without_values();
 
  // solution time series
  SmartPtr<vector_time_series_type> m_spSolTimeSeries;
  m_spSolTimeSeries=make_sp(new vector_time_series_type());
  m_spSolTimeSeries->clear();
  m_spSolTimeSeries->push(uold, t0);
 
  // create matrix operator
  SmartPtr<typename base_type::assembled_operator_type> spAssOp=make_sp(new typename base_type::assembled_operator_type(tdisc_dep_type::m_spTimeDisc, gl));
 
  // integrate
  if(!base_type::m_bNoLogOut)
    UG_LOG("+++ Integrating: ["<< t0 <<", "<< t1 <<"]\n");
 
  double t = t0;
  number dt_assembled = -1.0;   // invalid
  int step = 1;
 
  number currdt = base_type::m_dt;
 
  while((t < t1) && (t1-t > base_type::m_precisionBound))
  {
 
    if(!base_type::m_bNoLogOut)
    { UG_LOG("+++ Timestep +++" << step << "\n"); }
 
    // determine step size
    number dt = std::min(currdt, t1-t);
 
    // prepare step
    tdisc.prepare_step(m_spSolTimeSeries, dt);
    if (fabs(dt-dt_assembled) > base_type::m_precisionBound)
    {
      // re-assemble operator
      if(!base_type::m_bNoLogOut)
        UG_LOG("+++ Reassemble (t=" << t << ", dt=" << dt <<")\n");
 
      tdisc.assemble_linear(*spAssOp, *b, gl);
      (base_type::m_spLinearSolver)->init(spAssOp, *u1);
      dt_assembled = dt;
    }
    else
    {
      // keep old operator
      tdisc.assemble_rhs(*b, gl);
    }
 
    // execute step
    if (base_type::m_spLinearSolver->apply(*u1, *b))
    {
      // ACCEPTING:
      // Notify observers (which could modify u1!)
      this->notify_finalize_step(u1, step++, t+dt, dt);
 
      // prepare next step.
      // push updated solution into time series
      t += dt;
      SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries->oldest();
      VecAssign(*tmp, *u1.template cast_dynamic<typename base_type::vector_type>());
      m_spSolTimeSeries->push_discard_oldest(tmp, t);
 
    }
    else
    {
      // DISCARDING
      currdt *= 0.5;
    }
 
  }
 
  this->notify_end(u1, step, t1, currdt);
 
  return true;
};
 
 
template<typename TDomain, typename TAlgebra>
bool ConstStepLinearTimeIntegrator<TDomain, TAlgebra>::apply
(
  SmartPtr<grid_function_type> u1,
  number t1,
  ConstSmartPtr<grid_function_type> u0,
  number t0
)
{
  LIMEX_PROFILE_FUNC()
  UG_COND_THROW(!base_type::m_spLinearSolver.valid(), "Linear solver invalid");
 
  // short-cuts
  GridLevel const &gl = u0->grid_level();
  time_disc_type &tdisc = *tdisc_dep_type::m_spTimeDisc;
 
  // create solution vector & right hand side
  SmartPtr<typename base_type::vector_type> uold= u0->clone();
  SmartPtr<typename base_type::vector_type> b= u0->clone_without_values();
  
  // solution time series
  SmartPtr<vector_time_series_type> m_spSolTimeSeries;
  m_spSolTimeSeries=make_sp(new vector_time_series_type());
  m_spSolTimeSeries->push(uold, t0);
 
  SmartPtr<typename base_type::assembled_operator_type> spAssOp = SPNULL;
 
  // select number of steps
  double t = t0;
  int numSteps = round((t1-t0) / base_type::m_dt);
 
   // determine step size
  const double dt = (t1-t0) / numSteps;
  
   //std::cerr << "+++ Integrating: ["<< t0 <<", "<< t1 <<"] with dt=" << dt << "("<< numSteps<< " iters)\n";
  if(!base_type::m_bNoLogOut)
  {
    UG_LOG("+++ Integrating: [\t"<< t0 <<"\t, \t"<< t1 <<"\t] with dt=\t" << dt << "("<< numSteps<< " iters)" << std::endl);
  }
   
   // integrate
   for(int step = 1; step<=numSteps; ++step)
   {
 
    if(!base_type::m_bNoLogOut)
    {
      UG_LOG("+++ Const timestep +++" << step<< "(t=" << t << ", dt=" << dt << ")"<< std::endl);
    }
    this->notify_init_step(u1, step, t, dt);
    
    // prepare step
     tdisc.prepare_step(m_spSolTimeSeries, dt);
     if (spAssOp==SPNULL)
     {
       // Assemble operator.
      if(!base_type::m_bNoLogOut)
      { UG_LOG("+++ Assemble (t=" << t << ", dt=" << dt <<")" << std::endl); }
 
      spAssOp=make_sp(new typename base_type::assembled_operator_type(tdisc_dep_type::m_spTimeDisc, gl));
      tdisc.assemble_linear(*spAssOp, *b, gl);
      (base_type::m_spLinearSolver)->init(spAssOp, *u1);
     }
     else
     {
      // Recycle existing operator.
      // std::cerr << "Recycling timestep " << step << "\n";
      tdisc.assemble_rhs(*b, gl);
     }
 
     // execute step
     if (base_type::m_spLinearSolver->apply(*u1, *b))
     {
      // ACCEPTING:
      // Notify observers (which could modify u1!)
      this->notify_finalize_step(u1, step, t+dt, dt);
 
      // push updated solution into time series
      t += dt; // Update time.
      SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries->oldest();
      VecAssign(*tmp,  *u1.template cast_dynamic<typename base_type::vector_type>());
      m_spSolTimeSeries->push_discard_oldest(tmp, t);
      
     }
     else
     {
      UG_THROW("ConstStepLinearTimeIntegrator::apply failed!!!");
      this->notify_rewind_step(u1, step, t+dt, dt);
     }
 
   }
 
   this->notify_end(u1, numSteps, t1, dt);
 
   return true;
};
 
template<typename TDomain, typename TAlgebra>
bool TimeIntegratorLinearAdaptive<TDomain, TAlgebra>::apply
(
  SmartPtr<grid_function_type> u1,
  number t1,
  ConstSmartPtr<grid_function_type> u0,
  number t0
)
{
  // short-cuts
  GridLevel const &gl = u0->grid_level();
  time_disc_type &tdisc = *tdisc_dep_type::m_spTimeDisc;
  time_disc_type &tdisc2 = *m_spTimeDisc2;
 
  // create solution vector & right hand side
  SmartPtr<typename base_type::vector_type> uold= u0->clone();
  SmartPtr<typename base_type::vector_type> b= u0->clone_without_values();
 
  // create matrix operator
  SmartPtr<typename base_type::assembled_operator_type> spAssOp=make_sp(new typename base_type::assembled_operator_type(tdisc_dep_type::m_spTimeDisc, gl));
 
  // create additional solutions
  SmartPtr<typename base_type::grid_function_type>  u2old = u0->clone();
  SmartPtr<typename base_type::grid_function_type>  u2 = u0->clone();
 
  // solution time series
  SmartPtr<vector_time_series_type> m_spSolTimeSeries;
  m_spSolTimeSeries=make_sp(new vector_time_series_type());
  m_spSolTimeSeries->push(uold, t0);
 
  SmartPtr<vector_time_series_type> m_spSolTimeSeries2;
  m_spSolTimeSeries2=make_sp(new vector_time_series_type());
  m_spSolTimeSeries2->push(u2old, t0);
  
  // automatic selection of min/max time step
  if (m_dtmin <= 0.0) { m_dtmin = (t1 - t0)/1.0e+5; }
  if (m_dtmax <= 0.0) { m_dtmax = (t1 - t0)/10.0; }
 
  // Aitken Neville extrapolation
  const size_t tsteps[2] = {1,2};
  std::vector<size_t> vsteps (tsteps, tsteps+2);
  AitkenNevilleTimex<typename base_type::vector_type> timex(vsteps);
 
  // integrate
  if(!base_type::m_bNoLogOut)
    UG_LOG("+++ Integrating: ["<< t0 <<", "<< t1 <<"]\n");
 
   double t = t0;
   int step = 0;
 
   number dt = base_type::m_dt;
   while((t < t1) && (t1-t > base_type::m_precisionBound))
   {
     // step: t -> t+dt
     bool bSuccess = false;
     while (!bSuccess){
    
       // determine step size
       if (dt<m_dtmin){
        if(!base_type::m_bNoLogOut)
         UG_LOG("Step size below minimum")
     }
     dt = std::min(dt, t1-t);
     
     // basic step
      if(!base_type::m_bNoLogOut)
       UG_LOG("+++ Timestep: " << ++step << "\n");
 
     tdisc.prepare_step(m_spSolTimeSeries, dt);
     tdisc.assemble_linear(*spAssOp, *b, gl);
     base_type::m_spLinearSolver->init(spAssOp, *u1);
     base_type::m_spLinearSolver->apply(*u1, *b);
     
     
     // control 1/2
      if(!base_type::m_bNoLogOut)
       UG_LOG("+++ Control: " << step << "\n");
 
     tdisc2.prepare_step(m_spSolTimeSeries, 0.5*dt);
     tdisc2.assemble_linear(*spAssOp, *b, gl);
     base_type::m_spLinearSolver->init(spAssOp, *u2);
     base_type::m_spLinearSolver->apply(*u2, *b);
     
     // push back solution
     SmartPtr<typename base_type::vector_type> tmp2 =  m_spSolTimeSeries2->oldest();
     (*tmp2)=*(u2.template cast_static<typename base_type::vector_type>());
     m_spSolTimeSeries2->push_discard_oldest(tmp2, t + 0.5*dt);
     
     // control 2/2
     tdisc2.prepare_step(m_spSolTimeSeries2, 0.5*dt);
     tdisc2.assemble_linear(*spAssOp, *b, gl);
     base_type::m_spLinearSolver->init(spAssOp, *u2);
     base_type::m_spLinearSolver->apply(*u2, *b);
 
     // obtain extrapolated solution
     timex.set_solution(u1, 0);
     timex.set_solution(u2, 1);
     timex.apply();
 
     // predict (subsequent) time step
     number eps = timex.get_error_estimates()[0];
     number lambda = std::pow(0.8* m_tol/eps, 0.5); 
 
     number dtEst= dt*lambda; 
     dtEst = std::min(dtEst, 1.5*dt); 
     dtEst = std::min(dtEst, m_dtmax); 
     
     dt = dtEst;
     if (eps <= m_tol)
     {
       // ACCEPT STEP (and thus solution u2)
      if(!base_type::m_bNoLogOut)
       UG_LOG("ACCEPTING solution, dtnew=" << dt);
 
       bSuccess = true;
     }
     else
     {
       // DISCARD step
      if(!base_type::m_bNoLogOut)
       UG_LOG("DISCARDING solution, dtnew=" << dt);
 
       // => reset solutions
       VecAssign(*u1.template cast_dynamic<typename base_type::vector_type>(), *uold);
       
       // => timeSeries2 has been updated...
       SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries2->oldest();
       VecAssign(*tmp, *uold);
       m_spSolTimeSeries2->push_discard_oldest(tmp, t);
 
     }
 
     } 
 
   
     // prepare next loop
     t += dt;
 
     // push updated solution into time series (and continue)
     SmartPtr<typename base_type::vector_type> tmp = m_spSolTimeSeries->oldest();
     VecAssign(*tmp, static_cast<typename base_type::vector_type> (*u2));
     m_spSolTimeSeries->push_discard_oldest(tmp, t);
 
   }
 
   return true;
};
 
} // namespace ug
 
#endif /* __H__LIMEX__TIME_INTEGRATOR_IMPL_HPP__ */