Transient.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
 
#include "Transient.h"
 
// MOOSE includes
#include "Factory.h"
#include "SubProblem.h"
#include "TimeStepper.h"
#include "MooseApp.h"
#include "Conversion.h"
#include "FEProblem.h"
#include "NonlinearSystem.h"
#include "Control.h"
#include "TimePeriod.h"
#include "MooseMesh.h"
#include "AllLocalDofIndicesThread.h"
#include "TimeIntegrator.h"
#include "Console.h"
 
#include "libmesh/implicit_system.h"
#include "libmesh/nonlinear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/numeric_vector.h"
 
// C++ Includes
#include <iomanip>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
 
registerMooseObject("MooseApp", Transient);
 
defineLegacyParams(Transient);
 
InputParameters
Transient::validParams()
{
  InputParameters params = Executioner::validParams();
  params.addClassDescription("Executioner for time varying simulations.");
 
  std::vector<Real> sync_times(1);
  sync_times[0] = -std::numeric_limits<Real>::max();
 
  MooseEnum schemes("implicit-euler explicit-euler crank-nicolson bdf2 explicit-midpoint dirk "
                    "explicit-tvd-rk-2 newmark-beta",
                    "implicit-euler");
 
  params.addParam<Real>("start_time", 0.0, "The start time of the simulation");
  params.addParam<Real>("end_time", 1.0e30, "The end time of the simulation");
  params.addParam<Real>("dt", 1., "The timestep size between solves");
  params.addParam<Real>("dtmin", 2.0e-14, "The minimum timestep size in an adaptive run");
  params.addParam<Real>("dtmax", 1.0e30, "The maximum timestep size in an adaptive run");
  params.addParam<bool>(
      "reset_dt", false, "Use when restarting a calculation to force a change in dt.");
  params.addParam<unsigned int>("num_steps",
                                std::numeric_limits<unsigned int>::max(),
                                "The number of timesteps in a transient run");
  params.addParam<int>("n_startup_steps", 0, "The number of timesteps during startup");
 
  params.addDeprecatedParam<bool>("trans_ss_check",
                                  false,
                                  "Whether or not to check for steady state conditions",
                                  "Use steady_state_detection instead");
  params.addDeprecatedParam<Real>("ss_check_tol",
                                  1.0e-08,
                                  "Whenever the relative residual changes by less "
                                  "than this the solution will be considered to be "
                                  "at steady state.",
                                  "Use steady_state_tolerance instead");
  params.addDeprecatedParam<Real>(
      "ss_tmin",
      0.0,
      "Minimum amount of time to run before checking for steady state conditions.",
      "Use steady_state_start_time instead");
 
  params.addParam<bool>(
      "steady_state_detection", false, "Whether or not to check for steady state conditions");
  params.addParam<Real>("steady_state_tolerance",
                        1.0e-08,
                        "Whenever the relative residual changes by less "
                        "than this the solution will be considered to be "
                        "at steady state.");
  params.addParam<Real>(
      "steady_state_start_time",
      0.0,
      "Minimum amount of time to run before checking for steady state conditions.");
 
  params.addParam<std::vector<std::string>>("time_periods", "The names of periods");
  params.addParam<std::vector<Real>>("time_period_starts", "The start times of time periods");
  params.addParam<std::vector<Real>>("time_period_ends", "The end times of time periods");
  params.addParam<bool>(
      "abort_on_solve_fail", false, "abort if solve not converged rather than cut timestep");
  params.addParam<MooseEnum>("scheme", schemes, "Time integration scheme used.");
  params.addParam<Real>("timestep_tolerance",
                        2.0e-14,
                        "the tolerance setting for final timestep size and sync times");
 
  params.addParam<bool>("use_multiapp_dt",
                        false,
                        "If true then the dt for the simulation will be "
                        "chosen by the MultiApps.  If false (the "
                        "default) then the minimum over the master dt "
                        "and the MultiApps is used");
 
  params.addParamNamesToGroup(
      "steady_state_detection steady_state_tolerance steady_state_start_time",
      "Steady State Detection");
 
  params.addParamNamesToGroup("start_time dtmin dtmax n_startup_steps trans_ss_check ss_check_tol "
                              "ss_tmin abort_on_solve_fail timestep_tolerance use_multiapp_dt",
                              "Advanced");
 
  params.addParamNamesToGroup("time_periods time_period_starts time_period_ends", "Time Periods");
 
  return params;
}
 
Transient::Transient(const InputParameters & parameters)
  : Executioner(parameters),
    _problem(_fe_problem),
    _nl(_fe_problem.getNonlinearSystemBase()),
    _time_scheme(getParam<MooseEnum>("scheme").getEnum<Moose::TimeIntegratorType>()),
    _t_step(_problem.timeStep()),
    _time(_problem.time()),
    _time_old(_problem.timeOld()),
    _dt(_problem.dt()),
    _dt_old(_problem.dtOld()),
    _unconstrained_dt(declareRecoverableData<Real>("unconstrained_dt", -1)),
    _at_sync_point(declareRecoverableData<bool>("at_sync_point", false)),
    _last_solve_converged(declareRecoverableData<bool>("last_solve_converged", true)),
    _xfem_repeat_step(false),
    _end_time(getParam<Real>("end_time")),
    _dtmin(getParam<Real>("dtmin")),
    _dtmax(getParam<Real>("dtmax")),
    _num_steps(getParam<unsigned int>("num_steps")),
    _n_startup_steps(getParam<int>("n_startup_steps")),
    _steady_state_detection(getParam<bool>("steady_state_detection")),
    _steady_state_tolerance(getParam<Real>("steady_state_tolerance")),
    _steady_state_start_time(getParam<Real>("steady_state_start_time")),
    _sync_times(_app.getOutputWarehouse().getSyncTimes()),
    _abort(getParam<bool>("abort_on_solve_fail")),
    _time_interval(declareRecoverableData<bool>("time_interval", false)),
    _start_time(getParam<Real>("start_time")),
    _timestep_tolerance(getParam<Real>("timestep_tolerance")),
    _target_time(declareRecoverableData<Real>("target_time", -1)),
    _use_multiapp_dt(getParam<bool>("use_multiapp_dt")),
    _solution_change_norm(declareRecoverableData<Real>("solution_change_norm", 0.0)),
    _sln_diff(_nl.addVector("sln_diff", false, PARALLEL)),
    _final_timer(registerTimedSection("final", 1))
{
  _picard_solve.setInnerSolve(_feproblem_solve);
 
  // Handle deprecated parameters
  if (!parameters.isParamSetByAddParam("trans_ss_check"))
    _steady_state_detection = getParam<bool>("trans_ss_check");
 
  if (!parameters.isParamSetByAddParam("ss_check_tol"))
    _steady_state_tolerance = getParam<Real>("ss_check_tol");
 
  if (!parameters.isParamSetByAddParam("ss_tmin"))
    _steady_state_start_time = getParam<Real>("ss_tmin");
 
  _t_step = 0;
  _dt = 0;
  _next_interval_output_time = 0.0;
 
  // Either a start_time has been forced on us, or we want to tell the App about what our start time
  // is (in case anyone else is interested.
  if (_app.hasStartTime())
    _start_time = _app.getStartTime();
  else if (parameters.isParamSetByUser("start_time") && !_app.isRecovering())
    _app.setStartTime(_start_time);
 
  _time = _time_old = _start_time;
  _problem.transient(true);
 
  setupTimeIntegrator();
 
  if (_app.halfTransient()) // Cut timesteps and end_time in half...
  {
    _end_time /= 2.0;
    _num_steps /= 2.0;
 
    if (_num_steps == 0) // Always do one step in the first half
      _num_steps = 1;
  }
}
 
void
Transient::init()
{
  if (!_time_stepper.get())
  {
    InputParameters pars = _app.getFactory().getValidParams("ConstantDT");
    pars.set<SubProblem *>("_subproblem") = &_problem;
    pars.set<Transient *>("_executioner") = this;
 
    if (!_pars.isParamSetByAddParam("end_time") && !_pars.isParamSetByAddParam("num_steps") &&
        _pars.isParamSetByAddParam("dt"))
      pars.set<Real>("dt") = (getParam<Real>("end_time") - getParam<Real>("start_time")) /
                             static_cast<Real>(getParam<unsigned int>("num_steps"));
    else
      pars.set<Real>("dt") = getParam<Real>("dt");
 
    pars.set<bool>("reset_dt") = getParam<bool>("reset_dt");
    _time_stepper = _app.getFactory().create<TimeStepper>("ConstantDT", "TimeStepper", pars);
  }
 
  _problem.execute(EXEC_PRE_MULTIAPP_SETUP);
  _problem.initialSetup();
 
  if (_app.isRestarting())
  {
    if (_app.hasStartTime())
      _time = _time_old = _app.getStartTime();
    else
      _time_old = _time;
  }
 
  _time_stepper->init();
 
  if (_app.isRecovering()) // Recover case
  {
    if (_t_step == 0)
      mooseError("Internal error in Transient executioner: _t_step is equal to 0 while recovering "
                 "in init().");
 
    _dt_old = _dt;
  }
}
 
void
Transient::preExecute()
{
  _time_stepper->preExecute();
 
  if (!_app.isRecovering())
  {
    _t_step = 0;
    _dt = 0;
    _next_interval_output_time = 0.0;
    if (!_app.isRestarting())
      _time = _time_old = _start_time;
 
    _problem.outputStep(EXEC_INITIAL);
 
    computeDT();
    _dt = getDT();
    if (_dt == 0)
      mooseError("Time stepper computed zero time step size on initial which is not allowed.\n"
                 "1. If you are using an existing time stepper, double check the values in your "
                 "input file or report an error.\n"
                 "2. If you are developing a new time stepper, make sure that initial time step "
                 "size in your code is computed correctly.");
    _nl.getTimeIntegrator()->init();
  }
}
 
void
Transient::preStep()
{
  _time_stepper->preStep();
}
 
void
Transient::postStep()
{
  _time_stepper->postStep();
}
 
void
Transient::execute()
{
  preExecute();
 
  // Start time loop...
  while (keepGoing())
  {
    incrementStepOrReject();
    preStep();
    computeDT();
    takeStep();
    endStep();
    postStep();
  }
 
  if (lastSolveConverged())
  {
    _t_step++;
    if (_picard_solve.hasPicardIteration())
    {
      _problem.finishMultiAppStep(EXEC_TIMESTEP_BEGIN);
      _problem.finishMultiAppStep(EXEC_TIMESTEP_END);
    }
  }
 
  if (!_app.halfTransient())
  {
    TIME_SECTION(_final_timer);
    _problem.execMultiApps(EXEC_FINAL);
    _problem.finalizeMultiApps();
    _problem.execute(EXEC_FINAL);
    _problem.outputStep(EXEC_FINAL);
  }
 
  // This method is to finalize anything else we want to do on the problem side.
  _problem.postExecute();
 
  // This method can be overridden for user defined activities in the Executioner.
  postExecute();
}
 
void
Transient::computeDT()
{
  _time_stepper->computeStep(); // This is actually when DT gets computed
}
 
void
Transient::incrementStepOrReject()
{
  if (lastSolveConverged())
  {
    if (!_xfem_repeat_step)
    {
#ifdef LIBMESH_ENABLE_AMR
      if (_t_step != 0)
        _problem.adaptMesh();
#endif
 
      _time_old = _time;
      _t_step++;
 
      _problem.advanceState();
 
      if (_t_step == 1)
        return;
 
      /*
       * Call the multi-app executioners endStep and
       * postStep methods when doing Picard. We do not perform these calls for
       * loose coupling because Transient::endStep and Transient::postStep get
       * called from TransientMultiApp::solveStep in that case.
       */
      if (_picard_solve.hasPicardIteration())
      {
        _problem.finishMultiAppStep(EXEC_TIMESTEP_BEGIN);
        _problem.finishMultiAppStep(EXEC_TIMESTEP_END);
      }
 
      /*
       * Ensure that we increment the sub-application time steps so that
       * when dt selection is made in the master application, we are using
       * the correct time step information
       */
      _problem.incrementMultiAppTStep(EXEC_TIMESTEP_BEGIN);
      _problem.incrementMultiAppTStep(EXEC_TIMESTEP_END);
    }
  }
  else
  {
    _problem.restoreMultiApps(EXEC_TIMESTEP_BEGIN, true);
    _problem.restoreMultiApps(EXEC_TIMESTEP_END, true);
    _time_stepper->rejectStep();
    _time = _time_old;
  }
}
 
void
Transient::takeStep(Real input_dt)
{
  _dt_old = _dt;
 
  if (input_dt == -1.0)
    _dt = computeConstrainedDT();
  else
    _dt = input_dt;
 
  _time_stepper->preSolve();
 
  // Increment time
  _time = _time_old + _dt;
 
  _problem.timestepSetup();
 
  _problem.onTimestepBegin();
 
  for (MooseIndex(_num_grid_steps) step = 0; step <= _num_grid_steps; ++step)
  {
    _time_stepper->step();
    _xfem_repeat_step = _picard_solve.XFEMRepeatStep();
 
    _last_solve_converged = _time_stepper->converged();
 
    if (!lastSolveConverged())
    {
      _console << "Aborting as solve did not converge\n";
      break;
    }
 
    if (step != _num_grid_steps)
      _problem.uniformRefine();
  }
 
  if (!(_problem.haveXFEM() && _picard_solve.XFEMRepeatStep()))
  {
    if (lastSolveConverged())
      _time_stepper->acceptStep();
    else
      _time_stepper->rejectStep();
  }
 
  _time = _time_old;
 
  _time_stepper->postSolve();
 
  _solution_change_norm = relativeSolutionDifferenceNorm() / _dt;
 
  return;
}
 
void
Transient::endStep(Real input_time)
{
  if (input_time == -1.0)
    _time = _time_old + _dt;
  else
    _time = input_time;
 
  if (lastSolveConverged())
  {
    if (_xfem_repeat_step)
      _time = _time_old;
    else
    {
      _nl.getTimeIntegrator()->postStep();
 
      // Compute the Error Indicators and Markers
      _problem.computeIndicators();
      _problem.computeMarkers();
 
      // Perform the output of the current time step
      _problem.outputStep(EXEC_TIMESTEP_END);
 
      // output
      if (_time_interval && (_time + _timestep_tolerance >= _next_interval_output_time))
        _next_interval_output_time += _time_interval_output_interval;
    }
  }
}
 
Real
Transient::computeConstrainedDT()
{
  //  // If start up steps are needed
  //  if (_t_step == 1 && _n_startup_steps > 1)
  //    _dt = _input_dt/(double)(_n_startup_steps);
  //  else if (_t_step == 1+_n_startup_steps && _n_startup_steps > 1)
  //    _dt = _input_dt;
 
  Real dt_cur = _dt;
  std::ostringstream diag;
 
  // After startup steps, compute new dt
  if (_t_step > _n_startup_steps)
    dt_cur = getDT();
 
  else
  {
    diag << "Timestep < n_startup_steps, using old dt: " << std::setw(9) << std::setprecision(6)
         << std::setfill('0') << std::showpoint << std::left << _dt << " tstep: " << _t_step
         << " n_startup_steps: " << _n_startup_steps << std::endl;
  }
  _unconstrained_dt = dt_cur;
 
  if (_verbose)
    _console << diag.str();
 
  diag.str("");
  diag.clear();
 
  // Allow the time stepper to limit the time step
  _at_sync_point = _time_stepper->constrainStep(dt_cur);
 
  // Don't let time go beyond next time interval output if specified
  if ((_time_interval) && (_time + dt_cur + _timestep_tolerance >= _next_interval_output_time))
  {
    dt_cur = _next_interval_output_time - _time;
    _at_sync_point = true;
 
    diag << "Limiting dt for time interval output at time: " << std::setw(9) << std::setprecision(6)
         << std::setfill('0') << std::showpoint << std::left << _next_interval_output_time
         << " dt: " << std::setw(9) << std::setprecision(6) << std::setfill('0') << std::showpoint
         << std::left << dt_cur << std::endl;
  }
 
  // Adjust to a target time if set
  if (_target_time > 0 && _time + dt_cur + _timestep_tolerance >= _target_time)
  {
    dt_cur = _target_time - _time;
    _at_sync_point = true;
 
    diag << "Limiting dt for target time: " << std::setw(9) << std::setprecision(6)
         << std::setfill('0') << std::showpoint << std::left << _next_interval_output_time
         << " dt: " << std::setw(9) << std::setprecision(6) << std::setfill('0') << std::showpoint
         << std::left << dt_cur << std::endl;
  }
 
  // Constrain by what the multi apps are doing
  Real multi_app_dt = _problem.computeMultiAppsDT(EXEC_TIMESTEP_BEGIN);
  if (_use_multiapp_dt || multi_app_dt < dt_cur)
  {
    dt_cur = multi_app_dt;
    _at_sync_point = false;
    diag << "Limiting dt for MultiApps: " << std::setw(9) << std::setprecision(6)
         << std::setfill('0') << std::showpoint << std::left << dt_cur << std::endl;
  }
  multi_app_dt = _problem.computeMultiAppsDT(EXEC_TIMESTEP_END);
  if (multi_app_dt < dt_cur)
  {
    dt_cur = multi_app_dt;
    _at_sync_point = false;
    diag << "Limiting dt for MultiApps: " << std::setw(9) << std::setprecision(6)
         << std::setfill('0') << std::showpoint << std::left << dt_cur << std::endl;
  }
 
  if (_verbose)
    _console << diag.str();
 
  return dt_cur;
}
 
Real
Transient::getDT()
{
  return _time_stepper->getCurrentDT();
}
 
bool
Transient::keepGoing()
{
  bool keep_going = !_problem.isSolveTerminationRequested();
 
  // Check for stop condition based upon steady-state check flag:
  if (lastSolveConverged())
  {
    if (!_xfem_repeat_step)
    {
      if (_steady_state_detection == true && _time > _steady_state_start_time)
      {
        // Check solution difference relative norm against steady-state tolerance
        if (_solution_change_norm < _steady_state_tolerance)
        {
          _console << "Steady-State Solution Achieved at time: " << _time << std::endl;
          // Output last solve if not output previously by forcing it
          keep_going = false;
        }
        else // Keep going
        {
          // Print steady-state relative error norm
          _console << "Steady-State Relative Differential Norm: " << _solution_change_norm
                   << std::endl;
        }
      }
 
      // Check for stop condition based upon number of simulation steps and/or solution end time:
      if (static_cast<unsigned int>(_t_step) >= _num_steps)
        keep_going = false;
 
      if ((_time >= _end_time) || (fabs(_time - _end_time) <= _timestep_tolerance))
        keep_going = false;
    }
  }
  else if (_abort)
  {
    _console << "Aborting as solve did not converge and input selected to abort" << std::endl;
    keep_going = false;
  }
 
  return keep_going;
}
 
void
Transient::estimateTimeError()
{
}
 
bool
Transient::lastSolveConverged() const
{
  return _last_solve_converged;
}
 
void
Transient::postExecute()
{
  _time_stepper->postExecute();
}
 
void
Transient::setTargetTime(Real target_time)
{
  _target_time = target_time;
}
 
Real
Transient::getSolutionChangeNorm()
{
  return _solution_change_norm;
}
 
void
Transient::setupTimeIntegrator()
{
  if (_pars.isParamSetByUser("scheme") && _problem.hasTimeIntegrator())
    mooseError("You cannot specify time_scheme in the Executioner and independently add a "
               "TimeIntegrator to the system at the same time");
 
  if (!_problem.hasTimeIntegrator())
  {
    // backwards compatibility
    std::string ti_str;
    using namespace Moose;
 
    switch (_time_scheme)
    {
      case TI_IMPLICIT_EULER:
        ti_str = "ImplicitEuler";
        break;
      case TI_EXPLICIT_EULER:
        ti_str = "ExplicitEuler";
        break;
      case TI_CRANK_NICOLSON:
        ti_str = "CrankNicolson";
        break;
      case TI_BDF2:
        ti_str = "BDF2";
        break;
      case TI_EXPLICIT_MIDPOINT:
        ti_str = "ExplicitMidpoint";
        break;
      case TI_LSTABLE_DIRK2:
        ti_str = "LStableDirk2";
        break;
      case TI_EXPLICIT_TVD_RK_2:
        ti_str = "ExplicitTVDRK2";
        break;
      case TI_NEWMARK_BETA:
        ti_str = "NewmarkBeta";
        break;
      default:
        mooseError("Unknown scheme: ", _time_scheme);
        break;
    }
 
    InputParameters params = _app.getFactory().getValidParams(ti_str);
    _problem.addTimeIntegrator(ti_str, ti_str, params);
  }
}
 
std::string
Transient::getTimeStepperName()
{
  if (_time_stepper)
  {
    TimeStepper & ts = *_time_stepper;
    return demangle(typeid(ts).name());
  }
  else
    return std::string();
}
 
Real
Transient::relativeSolutionDifferenceNorm()
{
  const NumericVector<Number> & current_solution = *_nl.currentSolution();
  const NumericVector<Number> & old_solution = _nl.solutionOld();
 
  _sln_diff = current_solution;
  _sln_diff -= old_solution;
 
  return (_sln_diff.l2_norm() / current_solution.l2_norm());
}