www.mooseframework.org
DT2.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 // MOOSE includes
16 #include "DT2.h"
17 
18 #include "AuxiliarySystem.h"
19 #include "FEProblem.h"
20 #include "NonlinearSystemBase.h"
21 #include "TimeIntegrator.h"
22 
23 #include "libmesh/implicit_system.h"
24 #include "libmesh/nonlinear_implicit_system.h"
25 #include "libmesh/nonlinear_solver.h"
26 #include "libmesh/transient_system.h"
27 #include "libmesh/numeric_vector.h"
28 
29 // C++ Includes
30 #include <iomanip>
31 
32 template <>
35 {
37  params.addParam<Real>("dt", 1., "The initial time step size.");
38  params.addRequiredParam<Real>("e_tol", "Target error tolerance.");
39  params.addRequiredParam<Real>("e_max", "Maximum acceptable error.");
40  params.addParam<Real>("max_increase", 1.0e9, "Maximum ratio that the time step can increase.");
41 
42  return params;
43 }
44 
45 DT2::DT2(const InputParameters & parameters)
46  : TimeStepper(parameters),
47  _u_diff(NULL),
48  _u1(NULL),
49  _u2(NULL),
50  _u_saved(NULL),
51  _u_older_saved(NULL),
52  _aux1(NULL),
53  _aux_saved(NULL),
54  _aux_older_saved(NULL),
55  _error(0.),
56  _e_tol(getParam<Real>("e_tol")),
57  _e_max(getParam<Real>("e_max")),
58  _max_increase(getParam<Real>("max_increase"))
59 {
60 }
61 
62 void
64 {
66  System & nl_sys = _fe_problem.getNonlinearSystemBase().system();
67  _u1 = &nl_sys.add_vector("u1", true, GHOSTED);
68  _u2 = &nl_sys.add_vector("u2", false, GHOSTED);
69  _u_diff = &nl_sys.add_vector("u_diff", false, GHOSTED);
70 
71  _u_saved = &nl_sys.add_vector("u_saved", false, GHOSTED);
72  _u_older_saved = &nl_sys.add_vector("u_older_saved", false, GHOSTED);
73 
74  TransientExplicitSystem & aux_sys = _fe_problem.getAuxiliarySystem().sys();
75  _aux1 = &aux_sys.add_vector("aux1", true, GHOSTED);
76  _aux_saved = &aux_sys.add_vector("aux_saved", false, GHOSTED);
77  _aux_older_saved = &aux_sys.add_vector("aux_older_saved", false, GHOSTED);
78 }
79 
80 void
82 {
84  TransientExplicitSystem & aux_sys = _fe_problem.getAuxiliarySystem().sys();
85 
86  // save solution vectors
87  *_u_saved = *nl_sys.currentSolution();
88  *_u_older_saved = nl_sys.solutionOlder();
89 
90  *_aux_saved = *aux_sys.current_local_solution;
91  *_aux_older_saved = *aux_sys.older_local_solution;
92 
93  _u_saved->close();
94  _u_older_saved->close();
95  _aux_saved->close();
96  _aux_older_saved->close();
97 }
98 
99 void
101 {
103  System & nl_sys = nl.system();
104  TransientExplicitSystem & aux_sys = _fe_problem.getAuxiliarySystem().sys();
105 
106  // solve the problem with full dt
107  _fe_problem.solve();
108  _converged = nl.converged();
109  if (_converged)
110  {
111  // save the solution (for time step with dt)
112  *_u1 = *nl.currentSolution();
113  _u1->close();
114  *_aux1 = *aux_sys.current_local_solution;
115  _aux1->close();
116 
117  // take two steps with dt/2
118  _console << "Taking two dt/2 time steps" << std::endl;
119 
120  // restore solutions to the original state
121  *nl_sys.current_local_solution = *_u_saved;
122  *aux_sys.current_local_solution = *_aux_saved;
123  nl_sys.current_local_solution->close();
124  aux_sys.current_local_solution->close();
125 
126  _dt = getCurrentDT() / 2; // cut the time step in half
127  _time = _time_old + _dt;
128 
129  // 1. step
132 
133  _console << " - 1. step" << std::endl;
135  nl.solve();
136  _converged = nl.converged();
137 
138  if (_converged)
139  {
140  nl_sys.update();
141 
144 
145  _time += _dt;
146  // 2. step
149 
150  _console << " - 2. step" << std::endl;
152  nl.solve();
153  _converged = nl.converged();
154  if (_converged)
155  {
156  nl_sys.update();
157 
158  *_u2 = *nl_sys.current_local_solution;
159  _u2->close();
160 
161  // compute error
162  *_u_diff = *_u2;
163  *_u_diff -= *_u1;
164  _u_diff->close();
165 
166  _error = (_u_diff->l2_norm() / std::max(_u1->l2_norm(), _u2->l2_norm())) / getCurrentDT();
167 
168  // restore _dt to its original value
169  _dt = getCurrentDT();
170  }
171  else
172  {
173  // solve failed, restore _time
174  _time = _time_old;
175  }
176  }
177  else
178  {
179  // solve failed, restore _time
180  _time = _time_old;
181  }
182 
183  if (!_converged)
184  {
185  *nl_sys.current_local_solution = *_u1;
186  nl.solutionOld() = *_u1;
187  nl.solutionOlder() = *_u_saved;
188 
189  *aux_sys.current_local_solution = *_aux1;
190  *aux_sys.old_local_solution = *_aux1;
191  *aux_sys.older_local_solution = *_aux_saved;
192 
193  nl_sys.current_local_solution->close();
194  nl.solutionOld().close();
195  nl.solutionOlder().close();
196  aux_sys.current_local_solution->close();
197  aux_sys.old_local_solution->close();
198  aux_sys.older_local_solution->close();
199  }
200  }
201 }
202 
203 Real
205 {
206  return getParam<Real>("dt");
207 }
208 
209 Real
211 {
212  Real curr_dt = getCurrentDT();
213  Real new_dt =
214  curr_dt * std::pow(_e_tol / _error,
216  if (new_dt / curr_dt > _max_increase)
217  new_dt = curr_dt * _max_increase;
218 
219  return new_dt;
220 }
221 
222 void
224 {
225  if (_error >= _e_max)
226  _console << "DT2Transient: Marking last solve not converged since |U2-U1|/max(|U2|,|U1|) = "
227  << _error << " >= e_max." << std::endl;
228 
230  System & nl_sys = nl.system();
231  TransientExplicitSystem & aux_sys = _fe_problem.getAuxiliarySystem().sys();
232 
233  // recover initial state
234  *nl_sys.current_local_solution = *_u_saved;
235  nl.solutionOld() = *_u_saved;
237 
238  *aux_sys.current_local_solution = *_aux_saved;
239  *aux_sys.old_local_solution = *_aux_saved;
240  *aux_sys.older_local_solution = *_aux_older_saved;
241 
242  nl_sys.solution->close();
243  nl.solutionOld().close();
244  nl.solutionOlder().close();
245  aux_sys.solution->close();
246  aux_sys.old_local_solution->close();
247  aux_sys.older_local_solution->close();
248 }
249 
250 bool
252 {
253  if (!_converged)
254  return false;
255 
256  if (_error < _e_max)
257  return true;
258  else
259  return false;
260 }
NumericVector< Number > * _aux_older_saved
Definition: DT2.h:57
virtual const NumericVector< Number > *& currentSolution() override
The solution vector that is currently being operated on.
NumericVector< Number > * _aux_saved
Definition: DT2.h:57
NumericVector< Number > * _aux1
Definition: DT2.h:57
virtual bool converged()=0
Returns the convergence state.
Real _error
global relative time discretization error estimate
Definition: DT2.h:60
virtual void step() override
Take a time step.
Definition: DT2.C:100
TimeIntegrator * getTimeIntegrator()
NonlinearSystemBase & getNonlinearSystemBase()
virtual Real computeInitialDT() override
Called to compute _current_dt for the first timestep.
Definition: DT2.C:204
void setSolverDefaults(FEProblemBase &problem)
Definition: Moose.C:1233
virtual NumericVector< Number > & solutionOld()=0
Real & _time_old
Definition: TimeStepper.h:130
Base class for time stepping.
Definition: TimeStepper.h:31
virtual TransientExplicitSystem & sys()
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
NumericVector< Number > * _u_diff
Definition: DT2.h:55
Object is evaluated at the end of every time step.
Definition: MooseTypes.h:100
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:56
Real _max_increase
maximum increase ratio
Definition: DT2.h:66
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:124
void execute(const ExecFlagType &exec_type)
Convenience function for performing execution of MOOSE systems.
NumericVector< Number > * _u1
Definition: DT2.h:55
Real _e_tol
error tolerance
Definition: DT2.h:62
virtual int order()=0
InputParameters validParams< DT2 >()
Definition: DT2.C:34
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:144
DT2(const InputParameters &parameters)
Definition: DT2.C:45
InputParameters validParams< TimeStepper >()
Definition: TimeStepper.C:22
NumericVector< Number > * _u_saved
Definition: DT2.h:56
virtual void rejectStep() override
This gets called when time step is rejected.
Definition: DT2.C:223
virtual System & system() override
Get the reference to the libMesh system.
virtual void onTimestepBegin() override
Real _e_max
maximal error
Definition: DT2.h:64
virtual void preExecute()
Definition: TimeStepper.C:65
virtual void preExecute() override
Definition: DT2.C:63
virtual bool converged() override
If the time step converged.
Definition: DT2.C:251
virtual NumericVector< Number > & solutionOlder()=0
virtual Real computeDT() override
Called to compute _current_dt for a normal step.
Definition: DT2.C:210
Object is evaluated at the beginning of every time step.
Definition: MooseTypes.h:102
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:132
Real getCurrentDT()
Get the current_dt.
Definition: TimeStepper.h:89
Real & _time
Values from executioner.
Definition: TimeStepper.h:129
NumericVector< Number > * _u2
Definition: DT2.h:55
virtual void preSolve() override
Definition: DT2.C:81