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        3133 : PseudoTimestep::validParams()
      22             : {
      23        3133 :   InputParameters params = GeneralPostprocessor::validParams();
      24             : 
      25             :   /// Enum that is used to select the timestep-selection method.
      26       12532 :   MooseEnum method("SER RDM EXP", "SER");
      27             : 
      28       12532 :   method.addDocumentation("SER", "Switched evolution relaxation");
      29       12532 :   method.addDocumentation("EXP", "Exponential progression");
      30       12532 :   method.addDocumentation("RDM", "Residual difference method");
      31             : 
      32       12532 :   params.addRequiredParam<MooseEnum>(
      33             :       "method", method, "The method used for pseudotimestep timemarching");
      34       18798 :   params.addRequiredRangeCheckedParam<Real>("initial_dt", "initial_dt > 0", "Initial timestep");
      35       18798 :   params.addRequiredRangeCheckedParam<Real>(
      36             :       "alpha", "alpha > 0", "The parameter alpha used in the scaling of the timestep");
      37       12532 :   params.addParam<Real>("max_dt", "The largest timestep allowed");
      38        9399 :   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       12532 :   params.set<ExecFlagEnum>("execute_on") = {EXEC_TIMESTEP_END, EXEC_INITIAL};
      44        6266 :   params.suppressParameter<ExecFlagEnum>("execute_on");
      45             : 
      46        3133 :   params.addClassDescription("Computes pseudo-time steps for obtaining steady-state solutions "
      47             :                              "through a pseudo transient process.");
      48             : 
      49        6266 :   return params;
      50        6266 : }
      51             : 
      52          36 : PseudoTimestep::PseudoTimestep(const InputParameters & parameters)
      53             :   : GeneralPostprocessor(parameters),
      54          36 :     _method(getParam<MooseEnum>("method").getEnum<PseudotimeMethod>()),
      55          72 :     _initial_dt(getParam<Real>("initial_dt")),
      56          72 :     _alpha(getParam<Real>("alpha")),
      57          72 :     _bound(isParamValid("max_dt")),
      58          48 :     _max_dt(_bound ? getParam<Real>("max_dt") : std::numeric_limits<Real>::max()),
      59          72 :     _residual_norms_sequence(declareRestartableData<std::vector<Real>>("residual_norms_sequence")),
      60         108 :     _iterations_step_sequence(declareRestartableData<std::vector<Real>>("iterations_step_sequence"))
      61             : {
      62          36 :   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          36 :   if (_method == PseudotimeMethod::SER)
      67             :   {
      68          24 :     if (parameters.isParamValid("iterations_window"))
      69          36 :       _iterations_window = getParam<unsigned int>("iterations_window");
      70             :     else
      71           0 :       _iterations_window = 1;
      72             :   }
      73          48 :   else if (parameters.isParamValid("iterations_window"))
      74           0 :     mooseError("The iterations window can be only provided for the SER method.");
      75          36 : }
      76             : 
      77             : void
      78         223 : PseudoTimestep::initialize()
      79             : {
      80             :   // start with the max
      81         223 :   _dt = std::numeric_limits<Real>::max();
      82         223 : }
      83             : 
      84             : Real
      85         223 : PseudoTimestep::getValue() const
      86             : {
      87         223 :   return _dt;
      88             : }
      89             : 
      90             : Real
      91         190 : PseudoTimestep::currentResidualNorm() const
      92             : {
      93         190 :   NonlinearSystemBase & nl = _fe_problem.getNonlinearSystemBase(_sys.number());
      94         190 :   Real res_norm = 0.0;
      95         380 :   for (const auto var_num : make_range(nl.system().n_vars()))
      96             :   {
      97             :     auto var_res =
      98         190 :         nl.system().calculate_norm(nl.getResidualNonTimeVector(), var_num, libMesh::DISCRETE_L2);
      99         190 :     res_norm = res_norm + std::pow(var_res, 2);
     100             :   }
     101         190 :   res_norm = std::sqrt(res_norm);
     102         190 :   return res_norm;
     103             : }
     104             : 
     105             : void
     106         190 : PseudoTimestep::outputPseudoTimestep(Real curr_time) const
     107             : {
     108         190 :   unsigned int curr_step = _fe_problem.timeStep();
     109         190 :   Real step_dt = _fe_problem.dt();
     110             : 
     111         190 :   unsigned int res_size = _residual_norms_sequence.size();
     112             : 
     113         190 :   _console << "Current step " << curr_step << " and current dt " << step_dt
     114         190 :            << " for  steady state residual " << _residual_norms_sequence[res_size - 1]
     115         190 :            << " at time " << curr_time << std::endl;
     116         190 : }
     117             : 
     118             : Real
     119          45 : PseudoTimestep::timestepSER() const
     120             : {
     121          45 :   const unsigned int curr_step = _fe_problem.timeStep();
     122             : 
     123          45 :   const unsigned int steps_size = _iterations_step_sequence.size();
     124          45 :   const unsigned int res_size = _residual_norms_sequence.size();
     125          45 :   const unsigned int prev_steps = std::min(_iterations_window + 1, curr_step);
     126             : 
     127         135 :   const Real factor = std::pow(_residual_norms_sequence[res_size - prev_steps] /
     128          45 :                                    _residual_norms_sequence[res_size - 1],
     129          45 :                                _alpha);
     130             : 
     131          45 :   const Real update_dt = _iterations_step_sequence[steps_size - 1] * factor;
     132          45 :   return update_dt;
     133             : }
     134             : 
     135             : Real
     136          67 : PseudoTimestep::timestepRDM() const
     137             : {
     138          67 :   const unsigned int res_size = _residual_norms_sequence.size();
     139          67 :   const unsigned int steps_size = _iterations_step_sequence.size();
     140             : 
     141          67 :   Real exponent = 0.0;
     142             : 
     143          67 :   if (_residual_norms_sequence[res_size - 1] < _residual_norms_sequence[res_size - 2])
     144         134 :     exponent = (_residual_norms_sequence[res_size - 2] - _residual_norms_sequence[res_size - 1]) /
     145          67 :                _residual_norms_sequence[res_size - 2];
     146          67 :   const Real update_dt = _iterations_step_sequence[steps_size - 1] * std::pow(_alpha, exponent);
     147          67 :   return update_dt;
     148             : }
     149             : 
     150             : Real
     151          45 : PseudoTimestep::timestepEXP() const
     152             : {
     153          45 :   const Real factor = MathUtils::pow(_alpha, _fe_problem.timeStep() - 1);
     154             : 
     155          45 :   const Real update_dt = _initial_dt * factor;
     156          45 :   return update_dt;
     157             : }
     158             : 
     159             : void
     160         223 : PseudoTimestep::execute()
     161             : {
     162         223 :   TransientBase * transient = dynamic_cast<TransientBase *>(_app.getExecutioner());
     163             : 
     164             :   Real res_norm;
     165             :   Real curr_dt;
     166             :   Real update_dt;
     167             : 
     168         223 :   if (_current_execute_flag == EXEC_INITIAL)
     169          33 :     _dt = _initial_dt;
     170             : 
     171             :   // at the end of each timestep call the postprocessor to set values for dt
     172         223 :   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         190 :     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         190 :     _residual_norms_sequence.resize(curr_step);
     184         190 :     _iterations_step_sequence.resize(curr_step);
     185             : 
     186         190 :     res_norm = currentResidualNorm();
     187         190 :     _residual_norms_sequence[curr_step - 1] = res_norm;
     188             : 
     189         190 :     curr_dt = transient->getDT();
     190         190 :     _iterations_step_sequence[curr_step - 1] = curr_dt;
     191             : 
     192         190 :     _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         190 :     if (_fe_problem.timeStep() > 1)
     197             :     {
     198         157 :       update_dt = curr_dt;
     199         157 :       switch (_method)
     200             :       {
     201          45 :         case PseudotimeMethod::SER:
     202          45 :           update_dt = timestepSER();
     203          45 :           break;
     204          45 :         case PseudotimeMethod::EXP:
     205          45 :           update_dt = timestepEXP();
     206          45 :           break;
     207          67 :         case PseudotimeMethod::RDM:
     208          67 :           update_dt = timestepRDM();
     209          67 :           break;
     210             :       }
     211         157 :       if (_bound)
     212          45 :         _dt = std::min(_max_dt, update_dt);
     213             :       else
     214         112 :         _dt = update_dt;
     215             :     }
     216         190 :     outputPseudoTimestep(_fe_problem.time());
     217             :   }
     218         223 : }