www.mooseframework.org
DT2.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 // MOOSE includes
11 #include "DT2.h"
12 
13 #include "AuxiliarySystem.h"
14 #include "FEProblem.h"
15 #include "NonlinearSystemBase.h"
16 #include "TimeIntegrator.h"
17 
18 #include "libmesh/implicit_system.h"
19 #include "libmesh/nonlinear_implicit_system.h"
20 #include "libmesh/nonlinear_solver.h"
21 #include "libmesh/transient_system.h"
22 #include "libmesh/numeric_vector.h"
23 
24 // C++ Includes
25 #include <iomanip>
26 
27 registerMooseObject("MooseApp", DT2);
28 
29 template <>
32 {
34  params.addParam<Real>("dt", 1., "The initial time step size.");
35  params.addRequiredParam<Real>("e_tol", "Target error tolerance.");
36  params.addRequiredParam<Real>("e_max", "Maximum acceptable error.");
37  params.addParam<Real>("max_increase", 1.0e9, "Maximum ratio that the time step can increase.");
38 
39  return params;
40 }
41 
42 DT2::DT2(const InputParameters & parameters)
43  : TimeStepper(parameters),
44  _u_diff(NULL),
45  _u1(NULL),
46  _u2(NULL),
47  _u_saved(NULL),
48  _u_older_saved(NULL),
49  _aux1(NULL),
50  _aux_saved(NULL),
51  _aux_older_saved(NULL),
52  _error(0.),
53  _e_tol(getParam<Real>("e_tol")),
54  _e_max(getParam<Real>("e_max")),
55  _max_increase(getParam<Real>("max_increase"))
56 {
57 }
58 
59 void
61 {
63  System & nl_sys = _fe_problem.getNonlinearSystemBase().system();
64  _u1 = &nl_sys.add_vector("u1", true, GHOSTED);
65  _u2 = &nl_sys.add_vector("u2", false, GHOSTED);
66  _u_diff = &nl_sys.add_vector("u_diff", false, GHOSTED);
67 
68  _u_saved = &nl_sys.add_vector("u_saved", false, GHOSTED);
69  _u_older_saved = &nl_sys.add_vector("u_older_saved", false, GHOSTED);
70 
72  _aux1 = &aux_sys.add_vector("aux1", true, GHOSTED);
73  _aux_saved = &aux_sys.add_vector("aux_saved", false, GHOSTED);
74  _aux_older_saved = &aux_sys.add_vector("aux_older_saved", false, GHOSTED);
75 }
76 
77 void
79 {
82 
83  // save solution vectors
84  *_u_saved = *nl_sys.currentSolution();
85  *_u_older_saved = nl_sys.solutionOlder();
86 
87  *_aux_saved = *aux_sys.current_local_solution;
88  *_aux_older_saved = *aux_sys.older_local_solution;
89 
90  _u_saved->close();
91  _u_older_saved->close();
92  _aux_saved->close();
93  _aux_older_saved->close();
94 }
95 
96 void
98 {
100  System & nl_sys = nl.system();
102 
103  // solve the problem with full dt
104  _fe_problem.solve();
105  _converged = nl.converged();
106  if (_converged)
107  {
108  // save the solution (for time step with dt)
109  *_u1 = *nl.currentSolution();
110  _u1->close();
111  *_aux1 = *aux_sys.current_local_solution;
112  _aux1->close();
113 
114  // take two steps with dt/2
115  _console << "Taking two dt/2 time steps" << std::endl;
116 
117  // restore solutions to the original state
118  *nl_sys.current_local_solution = *_u_saved;
119  *aux_sys.current_local_solution = *_aux_saved;
120  nl_sys.current_local_solution->close();
121  aux_sys.current_local_solution->close();
122 
123  _dt = getCurrentDT() / 2; // cut the time step in half
124  _time = _time_old + _dt;
125 
126  // 1. step
129 
130  _console << " - 1. step" << std::endl;
132  nl.solve();
133  _converged = nl.converged();
134 
135  if (_converged)
136  {
137  nl_sys.update();
138 
141 
142  _time += _dt;
143  // 2. step
146 
147  _console << " - 2. step" << std::endl;
149  nl.solve();
150  _converged = nl.converged();
151  if (_converged)
152  {
153  nl_sys.update();
154 
155  *_u2 = *nl_sys.current_local_solution;
156  _u2->close();
157 
158  // compute error
159  *_u_diff = *_u2;
160  *_u_diff -= *_u1;
161  _u_diff->close();
162 
163  _error = (_u_diff->l2_norm() / std::max(_u1->l2_norm(), _u2->l2_norm())) / getCurrentDT();
164 
165  // restore _dt to its original value
166  _dt = getCurrentDT();
167  }
168  else
169  {
170  // solve failed, restore _time
171  _time = _time_old;
172  }
173  }
174  else
175  {
176  // solve failed, restore _time
177  _time = _time_old;
178  }
179 
180  if (!_converged)
181  {
182  *nl_sys.current_local_solution = *_u1;
183  nl.solutionOld() = *_u1;
185 
186  *aux_sys.current_local_solution = *_aux1;
187  *aux_sys.old_local_solution = *_aux1;
188  *aux_sys.older_local_solution = *_aux_saved;
189 
190  nl_sys.current_local_solution->close();
191  nl.solutionOld().close();
192  nl.solutionOlder().close();
193  aux_sys.current_local_solution->close();
194  aux_sys.old_local_solution->close();
195  aux_sys.older_local_solution->close();
196  }
197  }
198 }
199 
200 Real
202 {
203  return getParam<Real>("dt");
204 }
205 
206 Real
208 {
209  Real curr_dt = getCurrentDT();
210  Real new_dt =
211  curr_dt * std::pow(_e_tol / _error,
213  if (new_dt / curr_dt > _max_increase)
214  new_dt = curr_dt * _max_increase;
215 
216  return new_dt;
217 }
218 
219 void
221 {
222  if (_error >= _e_max)
223  _console << "DT2Transient: Marking last solve not converged since |U2-U1|/max(|U2|,|U1|) = "
224  << _error << " >= e_max." << std::endl;
225 
227  System & nl_sys = nl.system();
229 
230  // recover initial state
231  *nl_sys.current_local_solution = *_u_saved;
232  nl.solutionOld() = *_u_saved;
234 
235  *aux_sys.current_local_solution = *_aux_saved;
236  *aux_sys.old_local_solution = *_aux_saved;
237  *aux_sys.older_local_solution = *_aux_older_saved;
238 
239  nl_sys.solution->close();
240  nl.solutionOld().close();
241  nl.solutionOlder().close();
242  aux_sys.solution->close();
243  aux_sys.old_local_solution->close();
244  aux_sys.older_local_solution->close();
245 }
246 
247 bool
249 {
250  if (!_converged)
251  return false;
252 
253  if (_error < _e_max)
254  return true;
255  else
256  return false;
257 }
NumericVector< Number > * _aux_older_saved
Definition: DT2.h:51
NumericVector< Number > * _aux_saved
Definition: DT2.h:51
NumericVector< Number > * _aux1
Definition: DT2.h:51
virtual bool converged()=0
Returns the convergence state.
Real _error
global relative time discretization error estimate
Definition: DT2.h:54
virtual void step() override
Take a time step.
Definition: DT2.C:97
NonlinearSystemBase & getNonlinearSystemBase()
virtual Real computeInitialDT() override
Called to compute _current_dt for the first timestep.
Definition: DT2.C:201
void setSolverDefaults(FEProblemBase &problem)
Definition: Moose.C:537
virtual NumericVector< Number > & solutionOld()=0
Real & _time_old
Definition: TimeStepper.h:125
TimeIntegrator * getTimeIntegrator()
Definition: SystemBase.h:729
Base class for time stepping.
Definition: TimeStepper.h:26
virtual TransientExplicitSystem & sys()
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
virtual bool converged() const override
If the time step converged.
Definition: DT2.C:248
const ExecFlagType EXEC_TIMESTEP_END
NumericVector< Number > * _u_diff
Definition: DT2.h:49
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...
NumericVector< Number > * _u_older_saved
Definition: DT2.h:50
Real _max_increase
maximum increase ratio
Definition: DT2.h:60
An adaptive timestepper that compares the solution obtained from a single step of size dt with two st...
Definition: DT2.h:32
NonlinearSystemBase * nl
Nonlinear system to be solved.
virtual void advanceState()
Advance all of the state holding vectors / datastructures so that we can move to the next timestep...
FEProblemBase & _fe_problem
Definition: TimeStepper.h:119
virtual void execute(const ExecFlagType &exec_type)
Convenience function for performing execution of MOOSE systems.
NumericVector< Number > * _u1
Definition: DT2.h:49
Real pow(Real x, int e)
Definition: MathUtils.C:211
const ExecFlagType EXEC_TIMESTEP_BEGIN
registerMooseObject("MooseApp", DT2)
Real _e_tol
error tolerance
Definition: DT2.h:56
virtual int order()=0
TransientSystem< ExplicitSystem > TransientExplicitSystem
InputParameters validParams< DT2 >()
Definition: DT2.C:31
AuxiliarySystem & getAuxiliarySystem()
virtual void solve() override=0
Solve the system (using libMesh magic)
bool _converged
Whether or not the previous solve converged.
Definition: TimeStepper.h:139
DT2(const InputParameters &parameters)
Definition: DT2.C:42
const NumericVector< Number > *const & currentSolution() const override
The solution vector that is currently being operated on.
InputParameters validParams< TimeStepper >()
Definition: TimeStepper.C:17
NumericVector< Number > * _u_saved
Definition: DT2.h:50
virtual void rejectStep() override
This gets called when time step is rejected.
Definition: DT2.C:220
virtual System & system() override
Get the reference to the libMesh system.
virtual void onTimestepBegin() override
Real _e_max
maximal error
Definition: DT2.h:58
virtual void preExecute()
Definition: TimeStepper.C:61
virtual void preExecute() override
Definition: DT2.C:60
virtual NumericVector< Number > & solutionOlder()=0
virtual Real computeDT() override
Called to compute _current_dt for a normal step.
Definition: DT2.C:207
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...
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
virtual void solve() override
Real & _dt
Definition: TimeStepper.h:127
Real getCurrentDT()
Get the current_dt.
Definition: TimeStepper.h:84
Real & _time
Values from executioner.
Definition: TimeStepper.h:124
NumericVector< Number > * _u2
Definition: DT2.h:49
virtual void preSolve() override
Definition: DT2.C:78