IterationAdaptiveDT.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
// MOOSE includes
#include "IterationAdaptiveDT.h"
#include "Function.h"
#include "PiecewiseLinear.h"
#include "Transient.h"
#include "NonlinearSystem.h"
 
registerMooseObject("MooseApp", IterationAdaptiveDT);
 
defineLegacyParams(IterationAdaptiveDT);
 
InputParameters
IterationAdaptiveDT::validParams()
{
  InputParameters params = TimeStepper::validParams();
  params.addClassDescription("Adjust the timestep based on the number of iterations");
  params.addParam<int>("optimal_iterations",
                       "The target number of nonlinear iterations for adaptive timestepping");
  params.addParam<int>("iteration_window",
                       "Attempt to grow/shrink timestep if the iteration count "
                       "is below/above 'optimal_iterations plus/minus "
                       "iteration_window' (default = optimal_iterations/5).");
  params.addParam<unsigned>("linear_iteration_ratio",
                            "The ratio of linear to nonlinear iterations "
                            "to determine target linear iterations and "
                            "window for adaptive timestepping (default = "
                            "25)");
  params.addParam<PostprocessorName>("timestep_limiting_postprocessor",
                                     "If specified, the postprocessor value "
                                     "is used as an upper limit for the "
                                     "current time step length");
  params.addParam<FunctionName>("timestep_limiting_function",
                                "A 'PiecewiseBase' type function used to control the timestep by "
                                "limiting the change in the function over a timestep");
  params.addParam<Real>(
      "max_function_change",
      "The absolute value of the maximum change in timestep_limiting_function over a timestep");
  params.addParam<bool>("force_step_every_function_point",
                        false,
                        "Forces the timestepper to take "
                        "a step that is consistent with "
                        "points defined in the function");
  params.addRequiredParam<Real>("dt", "The default timestep size between solves");
  params.addParam<std::vector<Real>>("time_t", "The values of t");
  params.addParam<std::vector<Real>>("time_dt", "The values of dt");
  params.addParam<Real>("growth_factor",
                        2.0,
                        "Factor to apply to timestep if easy convergence (if "
                        "'optimal_iterations' is specified) or if recovering "
                        "from failed solve");
  params.addParam<Real>("cutback_factor",
                        0.5,
                        "Factor to apply to timestep if difficult "
                        "convergence (if 'optimal_iterations' is specified) "
                        "or if solution failed");
 
  params.addParam<bool>("reject_large_step",
                        false,
                        "If 'true', time steps that are too large compared to the "
                        "ideal time step will be rejected and repeated");
  params.addRangeCheckedParam<Real>("reject_large_step_threshold",
                                    0.1,
                                    "reject_large_step_threshold > 0 "
                                    "& reject_large_step_threshold < 1",
                                    "Ratio between the the ideal time step size and the "
                                    "current time step size below which a time step will "
                                    "be rejected if 'reject_large_step' is 'true'");
 
  params.declareControllable("growth_factor cutback_factor");
 
  return params;
}
 
IterationAdaptiveDT::IterationAdaptiveDT(const InputParameters & parameters)
  : TimeStepper(parameters),
    PostprocessorInterface(this),
    _dt_old(declareRestartableData<Real>("dt_old", 0.0)),
    _input_dt(getParam<Real>("dt")),
    _tfunc_last_step(declareRestartableData<bool>("tfunc_last_step", false)),
    _sync_last_step(declareRestartableData<bool>("sync_last_step", false)),
    _linear_iteration_ratio(isParamValid("linear_iteration_ratio")
                                ? getParam<unsigned>("linear_iteration_ratio")
                                : 25), // Default to 25
    _adaptive_timestepping(false),
    _pps_value(isParamValid("timestep_limiting_postprocessor")
                   ? &getPostprocessorValue("timestep_limiting_postprocessor")
                         : nullptr),
    _timestep_limiting_function(nullptr),
    _piecewise_timestep_limiting_function(nullptr),
    _piecewise_linear_timestep_limiting_function(nullptr),
    _times(0),
    _max_function_change(-1),
    _force_step_every_function_point(getParam<bool>("force_step_every_function_point")),
    _tfunc_times(getParam<std::vector<Real>>("time_t").begin(),
                 getParam<std::vector<Real>>("time_t").end()),
    _time_ipol(getParam<std::vector<Real>>("time_t"), getParam<std::vector<Real>>("time_dt")),
    _use_time_ipol(_time_ipol.getSampleSize() > 0),
    _growth_factor(getParam<Real>("growth_factor")),
    _cutback_factor(getParam<Real>("cutback_factor")),
    _nl_its(declareRestartableData<unsigned int>("nl_its", 0)),
    _l_its(declareRestartableData<unsigned int>("l_its", 0)),
    _cutback_occurred(declareRestartableData<bool>("cutback_occurred", false)),
    _at_function_point(false),
    _reject_large_step(getParam<bool>("reject_large_step")),
    _large_step_rejection_threshold(getParam<Real>("reject_large_step_threshold"))
{
  if (isParamValid("optimal_iterations"))
  {
    _adaptive_timestepping = true;
    _optimal_iterations = getParam<int>("optimal_iterations");
 
    if (isParamValid("iteration_window"))
      _iteration_window = getParam<int>("iteration_window");
    else
      _iteration_window = ceil(_optimal_iterations / 5.0);
  }
  else
  {
    if (isParamValid("iteration_window"))
      mooseError("'optimal_iterations' must be used for 'iteration_window' to be used");
    if (isParamValid("linear_iteration_ratio"))
      mooseError("'optimal_iterations' must be used for 'linear_iteration_ratio' to be used");
  }
 
  if (isParamValid("timestep_limiting_function"))
    _max_function_change =
        isParamValid("max_function_change") ? getParam<Real>("max_function_change") : -1;
  else
  {
    if (isParamValid("max_function_change"))
      mooseError("'timestep_limiting_function' must be used for 'max_function_change' to be used");
    if (_force_step_every_function_point)
      mooseError("'timestep_limiting_function' must be used for 'force_step_every_function_point' "
                 "to be used");
  }
}
 
void
IterationAdaptiveDT::init()
{
  if (isParamValid("timestep_limiting_function"))
  {
    _timestep_limiting_function =
        &_fe_problem.getFunction(getParam<FunctionName>("timestep_limiting_function"),
                                 isParamValid("_tid") ? getParam<THREAD_ID>("_tid") : 0);
    _piecewise_timestep_limiting_function =
        dynamic_cast<const PiecewiseBase *>(_timestep_limiting_function);
 
    if (_piecewise_timestep_limiting_function)
    {
      _piecewise_linear_timestep_limiting_function =
          dynamic_cast<const PiecewiseLinear *>(_piecewise_timestep_limiting_function);
      unsigned int time_size = _piecewise_timestep_limiting_function->functionSize();
      _times.resize(time_size);
 
      for (unsigned int i = 0; i < time_size; ++i)
        _times[i] = _piecewise_timestep_limiting_function->domain(i);
    }
    else
      mooseError("timestep_limiting_function must be a PiecewiseBase function");
  }
}
 
void
IterationAdaptiveDT::preExecute()
{
  TimeStepper::preExecute();
 
  // Delete all tfunc times that are at or before the begin time
  while (!_tfunc_times.empty() && _time + _timestep_tolerance >= *_tfunc_times.begin())
    _tfunc_times.erase(_tfunc_times.begin());
}
 
Real
IterationAdaptiveDT::computeInitialDT()
{
  return _input_dt;
}
 
Real
IterationAdaptiveDT::computeDT()
{
  Real dt = _dt_old;
 
  if (_cutback_occurred)
  {
    _cutback_occurred = false;
    if (_adaptive_timestepping)
    {
      // Don't allow it to grow this step, but shrink if needed
      bool allowToGrow = false;
      computeAdaptiveDT(dt, allowToGrow);
    }
  }
  else if (_tfunc_last_step)
  {
    _tfunc_last_step = false;
    _sync_last_step = false;
    dt = _time_ipol.sample(_time_old);
 
    if (_verbose)
    {
      _console << "Setting dt to value specified by dt function: " << std::setw(9) << dt << '\n';
    }
  }
  else if (_sync_last_step)
  {
    _sync_last_step = false;
    dt = _dt_old;
 
    if (_verbose)
    {
      _console << "Setting dt to value used before sync: " << std::setw(9) << dt << '\n';
    }
  }
  else if (_adaptive_timestepping)
    computeAdaptiveDT(dt);
  else if (_use_time_ipol)
    dt = computeInterpolationDT();
  else
  {
    dt *= _growth_factor;
    if (dt > _dt_old * _growth_factor)
      dt = _dt_old * _growth_factor;
  }
 
  return dt;
}
 
bool
IterationAdaptiveDT::constrainStep(Real & dt)
{
  bool at_sync_point = TimeStepper::constrainStep(dt);
 
  // Limit the timestep to postprocessor value
  limitDTToPostprocessorValue(dt);
 
  // Limit the timestep to limit change in the function
  limitDTByFunction(dt);
 
  // Adjust to the next tfunc time if needed
  if (!_tfunc_times.empty() && _time + dt + _timestep_tolerance >= *_tfunc_times.begin())
  {
    dt = *_tfunc_times.begin() - _time;
 
    if (_verbose)
    {
      _console << "Limiting dt to sync with dt function time: " << std::setw(9)
               << *_tfunc_times.begin() << " dt: " << std::setw(9) << dt << '\n';
    }
  }
 
  return at_sync_point;
}
 
Real
IterationAdaptiveDT::computeFailedDT()
{
  _cutback_occurred = true;
 
  // Can't cut back any more
  if (_dt <= _dt_min)
    mooseError("Solve failed and timestep already at dtmin, cannot continue!");
 
  if (_verbose)
  {
    _console << "\nSolve failed with dt: " << std::setw(9) << _dt
             << "\nRetrying with reduced dt: " << std::setw(9) << _dt * _cutback_factor << '\n';
  }
  else
    _console << "\nSolve failed, cutting timestep.\n";
 
  return _dt * _cutback_factor;
}
 
bool
IterationAdaptiveDT::converged() const
{
  if (!_reject_large_step)
    return TimeStepper::converged();
 
  // the solver has not converged
  if (!TimeStepper::converged())
    return false;
 
  // we are already at dt_min or at the start of the simulation
  // in which case we can move on to the next step
  if (_dt == _dt_min || _t_step < 2)
    return true;
 
  // we get what the next time step should be
  Real dt_test = _dt;
  limitDTToPostprocessorValue(dt_test);
 
  // we cannot constrain the time step any further
  if (dt_test == 0)
    return true;
 
  // if the time step is much smaller than the current time step
  // we need to repeat the current iteration with a smaller time step
  if (dt_test < _dt * _large_step_rejection_threshold)
    return false;
 
  // otherwise we move one
  return true;
}
 
void
IterationAdaptiveDT::limitDTToPostprocessorValue(Real & limitedDT) const
{
  if (_pps_value && _t_step > 1)
  {
    if (limitedDT > *_pps_value)
    {
      limitedDT = std::max(_dt_min, *_pps_value);
 
      if (_verbose)
        _console << "Limiting dt to postprocessor value. dt = " << limitedDT << '\n';
    }
  }
}
 
void
IterationAdaptiveDT::limitDTByFunction(Real & limitedDT)
{
  Real orig_dt = limitedDT;
 
  // Limit by function change for piecewise linear functions.
  if (_piecewise_linear_timestep_limiting_function && _max_function_change > 0)
  {
    const auto current_function_value =
        _piecewise_linear_timestep_limiting_function->value(_time, {});
 
    for (std::size_t next_time_index = 1; next_time_index < _times.size(); ++next_time_index)
    {
      const auto next_time = _times[next_time_index];
 
      // Skip ahead to find time point that is just past the current time.
      if (next_time + _timestep_tolerance <= _time)
        continue;
 
      // Find out how far we can go without exceeding the max function change.
      const auto next_function_value =
          _piecewise_linear_timestep_limiting_function->range(next_time_index);
      const auto change = std::abs(next_function_value - current_function_value);
      if (change > _max_function_change)
      {
        // Interpolate to find step.
        const auto restrictedStep = (_max_function_change / change) * (next_time - _time);
        if (restrictedStep < limitedDT)
        {
          limitedDT = std::max(_dt_min, restrictedStep);
        }
        break;
      }
 
      // Don't keep going if we've already passed the current limited step.
      if (next_time > _time + limitedDT)
        break;
    }
  }
  else if (_timestep_limiting_function)
  {
    Point dummyPoint;
    Real oldValue = _timestep_limiting_function->value(_time_old, dummyPoint);
    Real newValue = _timestep_limiting_function->value(_time_old + limitedDT, dummyPoint);
    Real change = std::abs(newValue - oldValue);
 
    if (_max_function_change > 0.0 && change > _max_function_change)
    {
      do
      {
        limitedDT /= 2.0;
        newValue = _timestep_limiting_function->value(_time_old + limitedDT, dummyPoint);
        change = std::abs(newValue - oldValue);
      } while (change > _max_function_change);
    }
  }
 
  _at_function_point = false;
  if (_piecewise_timestep_limiting_function && _force_step_every_function_point)
  {
    for (unsigned int i = 0; i + 1 < _times.size(); ++i)
    {
      if (_time >= _times[i] && _time < _times[i + 1])
      {
        if (limitedDT > _times[i + 1] - _time - _timestep_tolerance)
        {
          limitedDT = _times[i + 1] - _time;
          _at_function_point = true;
        }
        break;
      }
    }
  }
 
  if (_verbose && limitedDT != orig_dt)
  {
    if (_at_function_point)
      _console << "Limiting dt to match function point. dt = ";
    else
      _console << "Limiting dt to limit change in function. dt = ";
 
    _console << limitedDT << '\n';
  }
}
 
void
IterationAdaptiveDT::computeAdaptiveDT(Real & dt, bool allowToGrow, bool allowToShrink)
{
  const unsigned int growth_nl_its(
      _optimal_iterations > _iteration_window ? _optimal_iterations - _iteration_window : 0);
  const unsigned int shrink_nl_its(_optimal_iterations + _iteration_window);
  const unsigned int growth_l_its(_optimal_iterations > _iteration_window
                                      ? _linear_iteration_ratio *
                                            (_optimal_iterations - _iteration_window)
                                      : 0);
  const unsigned int shrink_l_its(_linear_iteration_ratio *
                                  (_optimal_iterations + _iteration_window));
 
  if (allowToGrow && (_nl_its < growth_nl_its && _l_its < growth_l_its))
  {
    // Grow the timestep
    dt *= _growth_factor;
 
    if (_verbose)
    {
      _console << "Growing dt: nl its = " << _nl_its << " < " << growth_nl_its
               << " && lin its = " << _l_its << " < " << growth_l_its << " old dt: " << std::setw(9)
               << _dt_old << " new dt: " << std::setw(9) << dt << '\n';
    }
  }
  else if (allowToShrink && (_nl_its > shrink_nl_its || _l_its > shrink_l_its))
  {
    // Shrink the timestep
    dt *= _cutback_factor;
 
    if (_verbose)
    {
      _console << "Shrinking dt: nl its = " << _nl_its << " > " << shrink_nl_its
               << " || lin its = " << _l_its << " > " << shrink_l_its << " old dt: " << std::setw(9)
               << _dt_old << " new dt: " << std::setw(9) << dt << '\n';
    }
  }
}
 
Real
IterationAdaptiveDT::computeInterpolationDT()
{
  Real dt = _time_ipol.sample(_time_old);
 
  if (dt > _dt_old * _growth_factor)
  {
    dt = _dt_old * _growth_factor;
 
    if (_verbose)
    {
      _console << "Growing dt to recover from cutback. "
               << " old dt: " << std::setw(9) << _dt_old << " new dt: " << std::setw(9) << dt
               << '\n';
    }
  }
 
  return dt;
}
 
void
IterationAdaptiveDT::rejectStep()
{
  TimeStepper::rejectStep();
}
 
void
IterationAdaptiveDT::acceptStep()
{
  TimeStepper::acceptStep();
 
  while (!_tfunc_times.empty() && _time + _timestep_tolerance >= *_tfunc_times.begin())
  {
    if (std::abs(_time - *_tfunc_times.begin()) <= _timestep_tolerance)
      _tfunc_last_step = true;
 
    _tfunc_times.erase(_tfunc_times.begin());
  }
 
  _nl_its = _fe_problem.getNonlinearSystemBase().nNonlinearIterations();
  _l_its = _fe_problem.getNonlinearSystemBase().nLinearIterations();
 
  if ((_at_function_point || _executioner.atSyncPoint()) &&
      _dt + _timestep_tolerance < _executioner.unconstrainedDT())
  {
    _dt_old = _fe_problem.dtOld();
    _sync_last_step = true;
 
    if (_verbose)
    {
      _console << "Sync point hit in current step, using previous dt for old dt: " << std::setw(9)
               << _dt_old << '\n';
    }
  }
  else
    _dt_old = _dt;
}