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