Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PseudoTimestep.h"
      11             : #include "FEProblemBase.h"
      12             : #include "NonlinearSystemBase.h"
      13             : #include "MathUtils.h"
      14             : #include "TransientBase.h"
      15             : #include "Restartable.h"
      16             : #include "libmesh/enum_norm_type.h"
      17             : 
      18             : registerMooseObject("MooseApp", PseudoTimestep);
      19             : 
      20             : InputParameters
      21       14343 : PseudoTimestep::validParams()
      22             : {
      23       14343 :   InputParameters params = GeneralPostprocessor::validParams();
      24             : 
      25             :   /// Enum that is used to select the timestep-selection method.
      26       14343 :   MooseEnum method("SER RDM EXP", "SER");
      27             : 
      28       14343 :   method.addDocumentation("SER", "Switched evolution relaxation");
      29       14343 :   method.addDocumentation("EXP", "Exponential progression");
      30       14343 :   method.addDocumentation("RDM", "Residual difference method");
      31             : 
      32       14343 :   params.addRequiredParam<MooseEnum>(
      33             :       "method", method, "The method used for pseudotimestep timemarching");
      34       14343 :   params.addRequiredRangeCheckedParam<Real>("initial_dt", "initial_dt > 0", "Initial timestep");
      35       14343 :   params.addRequiredRangeCheckedParam<Real>(
      36             :       "alpha", "alpha > 0", "The parameter alpha used in the scaling of the timestep");
      37       14343 :   params.addParam<Real>("max_dt", "The largest timestep allowed");
      38       14343 :   params.addParam<unsigned int>(
      39             :       "iterations_window",
      40             :       "For how many iterations should the residual be tracked (only applies to the SER method)");
      41             :   // Because this post-processor is meant to be used with PostprocessorDT, it
      42             :   // should be executed on initial (not included by default) and timestep end.
      43       43029 :   params.set<ExecFlagEnum>("execute_on") = {EXEC_TIMESTEP_END, EXEC_INITIAL};
      44       14343 :   params.suppressParameter<ExecFlagEnum>("execute_on");
      45             : 
      46       14343 :   params.addClassDescription("Computes pseudo-time steps for obtaining steady-state solutions "
      47             :                              "through a pseudo transient process.");
      48             : 
      49       28686 :   return params;
      50       28686 : }
      51             : 
      52          39 : PseudoTimestep::PseudoTimestep(const InputParameters & parameters)
      53             :   : GeneralPostprocessor(parameters),
      54          39 :     _method(getParam<MooseEnum>("method").getEnum<PseudotimeMethod>()),
      55          39 :     _initial_dt(getParam<Real>("initial_dt")),
      56          39 :     _alpha(getParam<Real>("alpha")),
      57          39 :     _bound(isParamValid("max_dt")),
      58          39 :     _max_dt(_bound ? getParam<Real>("max_dt") : std::numeric_limits<Real>::max()),
      59          39 :     _residual_norms_sequence(declareRestartableData<std::vector<Real>>("residual_norms_sequence")),
      60          78 :     _iterations_step_sequence(declareRestartableData<std::vector<Real>>("iterations_step_sequence"))
      61             : {
      62          39 :   if (!_fe_problem.isTransient())
      63           0 :     mooseError("This pseudotimestepper can only be used if the steady state is computed via time "
      64             :                "integration.");
      65             : 
      66          39 :   if (_method == PseudotimeMethod::SER)
      67             :   {
      68          13 :     if (parameters.isParamValid("iterations_window"))
      69          13 :       _iterations_window = getParam<unsigned int>("iterations_window");
      70             :     else
      71           0 :       _iterations_window = 1;
      72             :   }
      73          26 :   else if (parameters.isParamValid("iterations_window"))
      74           0 :     mooseError("The iterations window can be only provided for the SER method.");
      75          39 : }
      76             : 
      77             : void
      78         243 : PseudoTimestep::initialize()
      79             : {
      80             :   // start with the max
      81         243 :   _dt = std::numeric_limits<Real>::max();
      82         243 : }
      83             : 
      84             : Real
      85         243 : PseudoTimestep::getValue() const
      86             : {
      87         243 :   return _dt;
      88             : }
      89             : 
      90             : Real
      91         207 : PseudoTimestep::currentResidualNorm() const
      92             : {
      93         207 :   NonlinearSystemBase & nl = _fe_problem.getNonlinearSystemBase(_sys.number());
      94         207 :   Real res_norm = 0.0;
      95         414 :   for (const auto var_num : make_range(nl.system().n_vars()))
      96             :   {
      97             :     auto var_res =
      98         207 :         nl.system().calculate_norm(nl.getResidualNonTimeVector(), var_num, libMesh::DISCRETE_L2);
      99         207 :     res_norm = res_norm + std::pow(var_res, 2);
     100             :   }
     101         207 :   res_norm = std::sqrt(res_norm);
     102         207 :   return res_norm;
     103             : }
     104             : 
     105             : void
     106         207 : PseudoTimestep::outputPseudoTimestep(Real curr_time) const
     107             : {
     108         207 :   unsigned int curr_step = _fe_problem.timeStep();
     109         207 :   Real step_dt = _fe_problem.dt();
     110             : 
     111         207 :   unsigned int res_size = _residual_norms_sequence.size();
     112             : 
     113         207 :   _console << "Current step " << curr_step << " and current dt " << step_dt
     114         207 :            << " for  steady state residual " << _residual_norms_sequence[res_size - 1]
     115         207 :            << " at time " << curr_time << std::endl;
     116         207 : }
     117             : 
     118             : Real
     119          49 : PseudoTimestep::timestepSER() const
     120             : {
     121          49 :   const unsigned int curr_step = _fe_problem.timeStep();
     122             : 
     123          49 :   const unsigned int steps_size = _iterations_step_sequence.size();
     124          49 :   const unsigned int res_size = _residual_norms_sequence.size();
     125          49 :   const unsigned int prev_steps = std::min(_iterations_window + 1, curr_step);
     126             : 
     127         147 :   const Real factor = std::pow(_residual_norms_sequence[res_size - prev_steps] /
     128          49 :                                    _residual_norms_sequence[res_size - 1],
     129          49 :                                _alpha);
     130             : 
     131          49 :   const Real update_dt = _iterations_step_sequence[steps_size - 1] * factor;
     132          49 :   return update_dt;
     133             : }
     134             : 
     135             : Real
     136          73 : PseudoTimestep::timestepRDM() const
     137             : {
     138          73 :   const unsigned int res_size = _residual_norms_sequence.size();
     139          73 :   const unsigned int steps_size = _iterations_step_sequence.size();
     140             : 
     141          73 :   Real exponent = 0.0;
     142             : 
     143          73 :   if (_residual_norms_sequence[res_size - 1] < _residual_norms_sequence[res_size - 2])
     144         146 :     exponent = (_residual_norms_sequence[res_size - 2] - _residual_norms_sequence[res_size - 1]) /
     145          73 :                _residual_norms_sequence[res_size - 2];
     146          73 :   const Real update_dt = _iterations_step_sequence[steps_size - 1] * std::pow(_alpha, exponent);
     147          73 :   return update_dt;
     148             : }
     149             : 
     150             : Real
     151          49 : PseudoTimestep::timestepEXP() const
     152             : {
     153          49 :   const Real factor = MathUtils::pow(_alpha, _fe_problem.timeStep() - 1);
     154             : 
     155          49 :   const Real update_dt = _initial_dt * factor;
     156          49 :   return update_dt;
     157             : }
     158             : 
     159             : void
     160         243 : PseudoTimestep::execute()
     161             : {
     162         243 :   TransientBase * transient = dynamic_cast<TransientBase *>(_app.getExecutioner());
     163             : 
     164             :   Real res_norm;
     165             :   Real curr_dt;
     166             :   Real update_dt;
     167             : 
     168         243 :   if (_current_execute_flag == EXEC_INITIAL)
     169          36 :     _dt = _initial_dt;
     170             : 
     171             :   // at the end of each timestep call the postprocessor to set values for dt
     172         243 :   if (_current_execute_flag == EXEC_TIMESTEP_END)
     173             :   {
     174             :     // This is all in case a timestep fails and needs to be re-done
     175             :     // Otherwise this is a simply a push_back operation for the vectors
     176             :     mooseAssert(_fe_problem.timeStep() >= 1, "Should at least be on the first time step.");
     177         207 :     const std::size_t curr_step = _fe_problem.timeStep();
     178             :     mooseAssert(_residual_norms_sequence.size() <= curr_step &&
     179             :                     curr_step - _residual_norms_sequence.size() <= 1,
     180             :                 "Vector is improper length.");
     181             :     mooseAssert(_residual_norms_sequence.size() == _iterations_step_sequence.size(),
     182             :                 "Vectors should be same length.");
     183         207 :     _residual_norms_sequence.resize(curr_step);
     184         207 :     _iterations_step_sequence.resize(curr_step);
     185             : 
     186         207 :     res_norm = currentResidualNorm();
     187         207 :     _residual_norms_sequence[curr_step - 1] = res_norm;
     188             : 
     189         207 :     curr_dt = transient->getDT();
     190         207 :     _iterations_step_sequence[curr_step - 1] = curr_dt;
     191             : 
     192         207 :     _dt = curr_dt;
     193             : 
     194             :     // since some of the residual methods require previous residuals the pseudo timesteppers
     195             :     // start at the second step in the computation
     196         207 :     if (_fe_problem.timeStep() > 1)
     197             :     {
     198         171 :       update_dt = curr_dt;
     199         171 :       switch (_method)
     200             :       {
     201          49 :         case PseudotimeMethod::SER:
     202          49 :           update_dt = timestepSER();
     203          49 :           break;
     204          49 :         case PseudotimeMethod::EXP:
     205          49 :           update_dt = timestepEXP();
     206          49 :           break;
     207          73 :         case PseudotimeMethod::RDM:
     208          73 :           update_dt = timestepRDM();
     209          73 :           break;
     210             :       }
     211         171 :       if (_bound)
     212          49 :         _dt = std::min(_max_dt, update_dt);
     213             :       else
     214         122 :         _dt = update_dt;
     215             :     }
     216         207 :     outputPseudoTimestep(_fe_problem.time());
     217             :   }
     218         243 : }