www.mooseframework.org
AB2PredictorCorrector.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "AB2PredictorCorrector.h"
11 #include "AdamsPredictor.h"
12 #include "Problem.h"
13 #include "FEProblem.h"
14 #include "MooseApp.h"
15 #include "NonlinearSystem.h"
16 #include "AuxiliarySystem.h"
17 #include "TimeIntegrator.h"
18 #include "Conversion.h"
19 
20 #include "libmesh/nonlinear_solver.h"
21 #include "libmesh/numeric_vector.h"
22 
23 // C++ Includes
24 #include <iomanip>
25 #include <iostream>
26 #include <fstream>
27 
29 
32 {
34  params.addClassDescription(
35  "Implements second order Adams-Bashforth method for timestep calculation.");
36  params.addRequiredParam<Real>("e_tol", "Target error tolerance.");
37  params.addRequiredParam<Real>("e_max", "Maximum acceptable error.");
38  params.addRequiredParam<Real>("dt", "Initial time step size");
39  params.addParam<Real>("max_increase", 1.0e9, "Maximum ratio that the time step can increase.");
40  params.addParam<int>(
41  "steps_between_increase", 1, "the number of time steps before recalculating dt");
42  params.addParam<int>("start_adapting", 2, "when to start taking adaptive time steps");
43  params.addParam<Real>("scaling_parameter", .8, "scaling parameter for dt selection");
44  return params;
45 }
46 
48  : TimeStepper(parameters),
49  _u1(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).addVector("u1", true, GHOSTED)),
50  _aux1(_fe_problem.getAuxiliarySystem().addVector("aux1", true, GHOSTED)),
51  _pred1(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).addVector("pred1", true, GHOSTED)),
52  _dt_full(declareRestartableData<Real>("dt_full", 0)),
53  _error(declareRestartableData<Real>("error", 0)),
54  _e_tol(getParam<Real>("e_tol")),
55  _e_max(getParam<Real>("e_max")),
56  _max_increase(getParam<Real>("max_increase")),
57  _steps_between_increase(getParam<int>("steps_between_increase")),
58  _dt_steps_taken(declareRestartableData<int>("dt_steps_taken", 0)),
59  _start_adapting(getParam<int>("start_adapting")),
60  _my_dt_old(declareRestartableData<Real>("my_dt_old", 0)),
61  _infnorm(declareRestartableData<Real>("infnorm", 0)),
62  _scaling_parameter(getParam<Real>("scaling_parameter"))
63 {
64  Real predscale = 1.;
65  InputParameters params = _app.getFactory().getValidParams("AdamsPredictor");
66  params.set<Real>("scale") = predscale;
67  _fe_problem.addPredictor("AdamsPredictor", "adamspredictor", params);
68 }
69 
70 void
72 {
74 }
75 
76 void
78 {
79  // save dt
80  _dt_full = _dt;
81 }
82 
83 void
85 {
88 
89  _fe_problem.solve(/*nl_sys=*/0);
90  _converged = _fe_problem.converged(/*nl_sys=*/0);
91  if (_converged)
92  {
93  _u1 = *nl.currentSolution();
94  _u1.close();
95 
96  _aux1 = *aux.currentSolution();
97  _aux1.close();
98  if (_t_step >= _start_adapting)
99  {
100  // Calculate error if past the first solve
102 
104  _e_max = 1.1 * _e_tol * _infnorm;
105  _console << "Time Error Estimate: " << _error << std::endl;
106  }
107  else
108  {
109  // First time step is problematic, sure we converged but what does that mean? We don't know.
110  // Nor can we calculate the error on the first time step.
111  }
112  }
113 }
114 
115 bool
117 {
118  if (!_converged)
119  return false;
120  if (_error < _e_max)
121  return true;
122  else
123  return false;
124 }
125 
126 void
128 {
129  if (!converged())
130  _dt_steps_taken = 0;
131 
132  if (_error >= _e_max)
133  _console << "Marking last solve not converged " << _error << " " << _e_max << std::endl;
134 }
135 
136 Real
138 {
139  if (_t_step <= _start_adapting)
140  return _dt;
141 
142  _my_dt_old = _dt;
143 
144  _dt_steps_taken += 1;
146  {
147 
148  Real new_dt = _dt_full * _scaling_parameter * std::pow(_infnorm * _e_tol / _error, 1.0 / 3.0);
149 
150  if (new_dt / _dt_full > _max_increase)
151  new_dt = _dt_full * _max_increase;
152  _dt_steps_taken = 0;
153  return new_dt;
154  }
155 
156  return _dt;
157 }
158 
159 Real
161 {
162  return getParam<Real>("dt");
163 }
164 
165 Real
167 {
170  auto scheme = Moose::stringToEnum<Moose::TimeIntegratorType>(ti->type());
171  Real dt_old = _my_dt_old;
172  if (dt_old == 0)
173  dt_old = _dt;
174 
175  switch (scheme)
176  {
178  {
179  _pred1 *= -1;
180  _pred1 += solution;
181  Real calc = _dt * _dt * .5;
182  _pred1 *= calc;
183  return _pred1.l2_norm();
184  }
186  {
187  _pred1 -= solution;
188  _pred1 *= (_dt) / (3.0 * (_dt + dt_old));
189  return _pred1.l2_norm();
190  }
191  case Moose::TI_BDF2:
192  {
193  _pred1 *= -1.0;
194  _pred1 += solution;
195  Real topcalc = 2.0 * (_dt + dt_old) * (_dt + dt_old);
196  Real bottomcalc = 6.0 * _dt * _dt + 12.0 * _dt * dt_old + 5.0 * dt_old * dt_old;
197  _pred1 *= topcalc / bottomcalc;
198 
199  return _pred1.l2_norm();
200  }
201  default:
202  break;
203  }
204  return -1;
205 }
static InputParameters validParams()
Definition: TimeStepper.C:16
AB2PredictorCorrector(const InputParameters &parameters)
Real & _error
global relative time discretization error estimate
virtual Real computeInitialDT() override
Called to compute _current_dt for the first timestep.
int _steps_between_increase
steps to take before increasing dt
virtual bool converged(const unsigned int nl_sys_num)
Eventually we want to convert this virtual over to taking a nonlinear system number argument...
Definition: SubProblem.h:101
Real _e_tol
error tolerance
NumericVector< Number > & _aux1
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
TimeIntegrator * getTimeIntegrator()
Definition: SystemBase.h:869
InputParameters getValidParams(const std::string &name) const
Get valid parameters for the object.
Definition: Factory.C:67
Base class for time stepping.
Definition: TimeStepper.h:22
const NumericVector< Number > *const & currentSolution() const override
The solution vector that is currently being operated on.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
NumericVector< Number > & _pred1
virtual void preSolve() override
Real & _dt_full
dt of the big step
virtual void solve(const unsigned int nl_sys_num)
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
Factory & getFactory()
Retrieve a writable reference to the Factory associated with this App.
Definition: MooseApp.h:408
virtual void step() override
Take a time step.
Nonlinear system to be solved.
FEProblemBase & _fe_problem
Definition: TimeStepper.h:126
virtual void addPredictor(const std::string &type, const std::string &name, InputParameters &parameters)
virtual void postSolve() override
virtual Real l2_norm() const =0
virtual Real estimateTimeError(NumericVector< Number > &sol)
GHOSTED
const std::string & type() const
Get the type of this class.
Definition: MooseBase.h:50
Real _e_max
maximal error
NonlinearSystemBase & getNonlinearSystemBase(const unsigned int sys_num)
Real _max_increase
maximum increase ratio
MooseApp & _app
The MOOSE application this is associated with.
Definition: MooseBase.h:69
Real _scaling_parameter
scaling_parameter for time step selection, default is 0.8
AuxiliarySystem & getAuxiliarySystem()
registerMooseObject("MooseApp", AB2PredictorCorrector)
bool _converged
Whether or not the previous solve converged.
Definition: TimeStepper.h:146
virtual void close()=0
const NumericVector< Number > *const & currentSolution() const override
The solution vector that is currently being operated on.
NumericVector< Number > & _u1
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Base class for time integrators.
virtual bool converged() const override
If the time step converged.
virtual NumericVector< Number > & solutionPredictor()
Definition: Predictor.h:41
virtual void preExecute()
Definition: TimeStepper.C:70
Real & _infnorm
infinity norm of the solution vector
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
int & _t_step
Definition: TimeStepper.h:133
virtual Real computeDT() override
Called to compute _current_dt for a normal step.
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
virtual void preExecute() override
int & _dt_steps_taken
steps taken at current dt
A TimeStepper based on the AB2 method.
Real & _dt
Definition: TimeStepper.h:134
MooseUnits pow(const MooseUnits &, int)
Definition: Units.C:537
A system that holds auxiliary variables.
virtual Real linfty_norm() const =0
void ErrorVector unsigned int
static InputParameters validParams()