Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_Stepper_impl.hpp
Go to the documentation of this file.
1 //@HEADER
2 // *****************************************************************************
3 // Tempus: Time Integration and Sensitivity Analysis Package
4 //
5 // Copyright 2017 NTESS and the Tempus contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 //@HEADER
9 
10 #ifndef Tempus_Stepper_impl_hpp
11 #define Tempus_Stepper_impl_hpp
12 
13 #include "NOX_Thyra.H"
14 
15 namespace Tempus {
16 
17 template <class Scalar>
19 {
20  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
21  out->setOutputToRootOnly(0);
22 
23  bool isValidSetup = this->isValidSetup(*out);
24 
25  if (isValidSetup)
26  this->isInitialized_ = true; // Only place it is set to true.
27  else
28  this->describe(*out, Teuchos::VERB_MEDIUM);
29 }
30 
31 template <class Scalar>
33 {
34  if (!this->isInitialized()) {
35  this->describe(*(this->getOStream()), Teuchos::VERB_MEDIUM);
37  !this->isInitialized(), std::logic_error,
38  "Error - " << this->description() << " is not initialized!");
39  }
40 }
41 
42 template <class Scalar>
44 {
45  if (a == false) {
46  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
47  out->setOutputToRootOnly(0);
48  Teuchos::OSTab ostab(out, 1, "Stepper::setUseFSALTrueOnly()");
49  *out << "Warning -- useFSAL for '" << this->getStepperType() << "'\n"
50  << "can only be set to true. Leaving set to true." << std::endl;
51  }
52  useFSAL_ = true;
53 }
54 
55 template <class Scalar>
57 {
58  if (a == true) {
59  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
60  out->setOutputToRootOnly(0);
61  Teuchos::OSTab ostab(out, 1, "Stepper::setUseFSALFalseOnly()");
62  *out << "Warning -- useFSAL for '" << this->getStepperType() << "'\n"
63  << "can only be set to false. Leaving set to false." << std::endl;
64  }
65  useFSAL_ = false;
66 }
67 
68 template <class Scalar>
70 {
72  stepperX_ == Teuchos::null, std::logic_error,
73  "Error - stepperX_ has not been set in setInitialConditions() or\n"
74  " can not be set from the state!\n");
75 
76  return stepperX_;
77 }
78 
79 template <class Scalar>
81 {
83  stepperXDot_ == Teuchos::null, std::logic_error,
84  "Error - stepperXDot_ has not been set in setInitialConditions() or\n"
85  " can not be set from the state!\n");
86 
87  return stepperXDot_;
88 }
89 
90 template <class Scalar>
92 {
94  stepperXDotDot_ == Teuchos::null, std::logic_error,
95  "Error - stepperXDotDot_ has not been set in setInitialConditions() or\n"
96  " can not be set from the state!\n");
97 
98  return stepperXDotDot_;
99 }
100 
101 template <class Scalar>
104 {
105  if (state->getXDotDot() != Teuchos::null)
106  stepperXDotDot_ = state->getXDotDot();
107  // Else use temporary storage stepperXDotDot_ which should have been set in
108  // setInitialConditions().
109 
111  stepperXDotDot_ == Teuchos::null, std::logic_error,
112  "Error - stepperXDotDot_ has not been set in setInitialConditions() or\n"
113  " can not be set from the state!\n");
114 
115  return stepperXDotDot_;
116 }
117 
118 template <class Scalar>
120  const Teuchos::EVerbosityLevel verbLevel) const
121 {
122  auto l_out = Teuchos::fancyOStream(out.getOStream());
123  Teuchos::OSTab ostab(*l_out, 2, this->description());
124  l_out->setOutputToRootOnly(0);
125 
126  *l_out << "--- Stepper ---\n"
127  << " isInitialized_ = " << Teuchos::toString(isInitialized_)
128  << std::endl
129  << " stepperType_ = " << stepperType_ << std::endl
130  << " useFSAL_ = " << Teuchos::toString(useFSAL_)
131  << std::endl
132  << " ICConsistency_ = " << ICConsistency_ << std::endl
133  << " ICConsistencyCheck_ = " << Teuchos::toString(ICConsistencyCheck_)
134  << std::endl
135  << " stepperX_ = " << stepperX_ << std::endl
136  << " stepperXDot_ = " << stepperXDot_ << std::endl
137  << " stepperXDotDot_ = " << stepperXDotDot_ << std::endl;
138 }
139 
140 template <class Scalar>
142 {
143  out.setOutputToRootOnly(0);
144  bool isValidSetup = true;
145 
146  if (!(ICConsistency_ == "None" || ICConsistency_ == "Zero" ||
147  ICConsistency_ == "App" || ICConsistency_ == "Consistent")) {
148  isValidSetup = false;
149  auto l_out = Teuchos::fancyOStream(out.getOStream());
150  l_out->setOutputToRootOnly(0);
151  *l_out << "The IC consistency does not have a valid value!\n"
152  << "('None', 'Zero', 'App' or 'Consistent')\n"
153  << " ICConsistency = " << ICConsistency_ << "\n";
154  }
155 
156  return isValidSetup;
157 }
158 
159 template <class Scalar>
161 {
162  if (pl != Teuchos::null) {
163  this->setStepperName(pl->name());
164  auto stepperType =
165  pl->get<std::string>("Stepper Type", this->getStepperType());
167  stepperType != this->getStepperType(), std::runtime_error,
168  " ParameterList 'Stepper Type' (='"
169  << stepperType
170  << "')\n does not match type for this Stepper (='"
171  << this->getStepperType() << "').");
172  this->setStepperType(stepperType);
173 
174  this->setUseFSAL(pl->get<bool>("Use FSAL", this->getUseFSAL()));
175  this->setICConsistency(pl->get<std::string>("Initial Condition Consistency",
176  this->getICConsistency()));
177  this->setICConsistencyCheck(pl->get<bool>(
178  "Initial Condition Consistency Check", this->getICConsistencyCheck()));
179  }
180 }
181 
182 template <class Scalar>
184  const
185 {
186  return this->getValidParametersBasic();
187 }
188 
189 template <class Scalar>
191  const
192 {
193  auto pl = Teuchos::parameterList(this->getStepperName());
194  pl->template set<std::string>("Stepper Type", this->getStepperType());
195 
196  pl->template set<bool>(
197  "Use FSAL", this->getUseFSAL(),
198  "The First-Same-As-Last (FSAL) principle is the situation where the\n"
199  "last function evaluation, f(x^{n-1},t^{n-1}) [a.k.a. xDot^{n-1}],\n"
200  "can be used for the first function evaluation, f(x^n,t^n)\n"
201  "[a.k.a. xDot^n]. For RK methods, this applies to the stages.\n"
202  "\n"
203  "Often the FSAL priniciple can be used to save an evaluation.\n"
204  "However there are cases when it cannot be used, e.g., operator\n"
205  "splitting where other steppers/operators have modified the solution,\n"
206  "x^*, and thus require the function evaluation, f(x^*, t^{n-1}).\n"
207  "\n"
208  "It should be noted that when the FSAL priniciple can be used\n"
209  "(can set useFSAL=true), setting useFSAL=false will give the\n"
210  "same solution but at additional expense. However, the reverse\n"
211  "is not true. When the FSAL priniciple can not be used\n"
212  "(need to set useFSAL=false), setting useFSAL=true will produce\n"
213  "incorrect solutions.\n"
214  "\n"
215  "Default in general for explicit and implicit steppers is false,\n"
216  "but individual steppers can override this default.");
217 
218  pl->template set<std::string>(
219  "Initial Condition Consistency", this->getICConsistency(),
220  "This indicates which type of consistency should be applied to\n"
221  "the initial conditions (ICs):\n"
222  "\n"
223  " 'None' - Do nothing to the ICs provided in the SolutionHistory.\n"
224  " 'Zero' - Set the derivative of the SolutionState to zero in the\n"
225  " SolutionHistory provided, e.g., xDot^0 = 0, or \n"
226  " xDotDot^0 = 0.\n"
227  " 'App' - Use the application's ICs, e.g., getNominalValues().\n"
228  " 'Consistent' - Make the initial conditions for x and xDot\n"
229  " consistent with the governing equations, e.g.,\n"
230  " xDot = f(x,t), and f(x, xDot, t) = 0. For implicit\n"
231  " ODEs, this requires a solve of f(x, xDot, t) = 0 for\n"
232  " xDot, and another Jacobian and residual may be\n"
233  " needed, e.g., boundary conditions on xDot may need\n"
234  " to replace boundary conditions on x.\n"
235  "\n"
236  "In general for explicit steppers, the default is 'Consistent',\n"
237  "because it is fairly cheap with just one residual evaluation.\n"
238  "In general for implicit steppers, the default is 'None', because\n"
239  "the application often knows its IC and can set it the initial\n"
240  "SolutionState. Also, as noted above, 'Consistent' may require\n"
241  "another Jacobian from the application. Individual steppers may\n"
242  "override these defaults.");
243 
244  pl->template set<bool>(
245  "Initial Condition Consistency Check", this->getICConsistencyCheck(),
246  "Check if the initial condition, x and xDot, is consistent with the\n"
247  "governing equations, xDot = f(x,t), or f(x, xDot, t) = 0.\n"
248  "\n"
249  "In general for explicit and implicit steppers, the default is true,\n"
250  "because it is fairly cheap with just one residual evaluation.\n"
251  "Individual steppers may override this default.");
252 
253  return pl;
254 }
255 
256 // Nonmember Helper Functions.
257 // ------------------------------------------------------------------------
258 
259 template <class Scalar>
261  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model)
262 {
264  typedef Thyra::ModelEvaluatorBase MEB;
265  const MEB::InArgs<Scalar> inArgs = model->createInArgs();
266  const MEB::OutArgs<Scalar> outArgs = model->createOutArgs();
267  const bool supports =
268  inArgs.supports(MEB::IN_ARG_x) && outArgs.supports(MEB::OUT_ARG_f);
269 
271  supports == false, std::logic_error,
272  model->description()
273  << " can not support an explicit ODE with\n"
274  << " IN_ARG_x = " << inArgs.supports(MEB::IN_ARG_x) << "\n"
275  << " OUT_ARG_f = " << outArgs.supports(MEB::OUT_ARG_f) << "\n"
276  << "Explicit ODE requires:\n"
277  << " IN_ARG_x = true\n"
278  << " OUT_ARG_f = true\n"
279  << "\n"
280  << "NOTE: Currently the convention to evaluate f(x,t) is to set\n"
281  << "xdot=null! There is no InArgs support to test if xdot is null,\n"
282  << "so we set xdot=null and hope the ModelEvaluator can handle "
283  "it.\n");
284 
285  return;
286 }
287 
288 template <class Scalar>
290  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model)
291 {
293  typedef Thyra::ModelEvaluatorBase MEB;
294  const MEB::InArgs<Scalar> inArgs = model->createInArgs();
295  const MEB::OutArgs<Scalar> outArgs = model->createOutArgs();
296  const bool supports = inArgs.supports(MEB::IN_ARG_x) &&
297  inArgs.supports(MEB::IN_ARG_x_dot) &&
298  outArgs.supports(MEB::OUT_ARG_f);
299 
301  supports == false, std::logic_error,
302  model->description()
303  << "can not support an explicit ODE with\n"
304  << " IN_ARG_x = " << inArgs.supports(MEB::IN_ARG_x) << "\n"
305  << " IN_ARG_x_dot = " << inArgs.supports(MEB::IN_ARG_x_dot) << "\n"
306  << " OUT_ARG_f = " << outArgs.supports(MEB::OUT_ARG_f) << "\n"
307  << "Explicit ODE requires:\n"
308  << " IN_ARG_x = true\n"
309  << " IN_ARG_x_dot = true\n"
310  << " OUT_ARG_f = true\n"
311  << "\n"
312  << "NOTE: Currently the convention to evaluate f(x, xdot, t) is to\n"
313  << "set xdotdot=null! There is no InArgs support to test if "
314  "xdotdot\n"
315  << "is null, so we set xdotdot=null and hope the ModelEvaluator can\n"
316  << "handle it.\n");
317 
318  return;
319 }
320 
321 template <class Scalar>
323  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model)
324 {
326  typedef Thyra::ModelEvaluatorBase MEB;
327  const MEB::InArgs<Scalar> inArgs = model->createInArgs();
328  const MEB::OutArgs<Scalar> outArgs = model->createOutArgs();
329  const bool supports =
330  inArgs.supports(MEB::IN_ARG_x) && inArgs.supports(MEB::IN_ARG_x_dot) &&
331  inArgs.supports(MEB::IN_ARG_alpha) && inArgs.supports(MEB::IN_ARG_beta) &&
332  !inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff) &&
333  outArgs.supports(MEB::OUT_ARG_f) && outArgs.supports(MEB::OUT_ARG_W);
334 
336  supports == false, std::logic_error,
337  model->description() << " can not support an implicit ODE with\n"
338  << " IN_ARG_x = "
339  << inArgs.supports(MEB::IN_ARG_x) << "\n"
340  << " IN_ARG_x_dot = "
341  << inArgs.supports(MEB::IN_ARG_x_dot) << "\n"
342  << " IN_ARG_alpha = "
343  << inArgs.supports(MEB::IN_ARG_alpha) << "\n"
344  << " IN_ARG_beta = "
345  << inArgs.supports(MEB::IN_ARG_beta) << "\n"
346  << " IN_ARG_W_x_dot_dot_coeff = "
347  << inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)
348  << "\n"
349  << " OUT_ARG_f = "
350  << outArgs.supports(MEB::OUT_ARG_f) << "\n"
351  << " OUT_ARG_W = "
352  << outArgs.supports(MEB::OUT_ARG_W) << "\n"
353  << "Implicit ODE requires:\n"
354  << " IN_ARG_x = true\n"
355  << " IN_ARG_x_dot = true\n"
356  << " IN_ARG_alpha = true\n"
357  << " IN_ARG_beta = true\n"
358  << " IN_ARG_W_x_dot_dot_coeff = false\n"
359  << " OUT_ARG_f = true\n"
360  << " OUT_ARG_W = true\n");
361 
362  return;
363 }
364 
365 template <class Scalar>
367  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model)
368 {
370  typedef Thyra::ModelEvaluatorBase MEB;
371  const MEB::InArgs<Scalar> inArgs = model->createInArgs();
372  const MEB::OutArgs<Scalar> outArgs = model->createOutArgs();
373  const bool supports =
374  inArgs.supports(MEB::IN_ARG_x) && inArgs.supports(MEB::IN_ARG_x_dot) &&
375  inArgs.supports(MEB::IN_ARG_x_dot_dot) &&
376  inArgs.supports(MEB::IN_ARG_alpha) && inArgs.supports(MEB::IN_ARG_beta) &&
377  inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff) &&
378  outArgs.supports(MEB::OUT_ARG_f) && outArgs.supports(MEB::OUT_ARG_W);
379 
381  supports == false, std::logic_error,
382  model->description() << " can not support an implicit ODE with\n"
383  << " IN_ARG_x = "
384  << inArgs.supports(MEB::IN_ARG_x) << "\n"
385  << " IN_ARG_x_dot = "
386  << inArgs.supports(MEB::IN_ARG_x_dot) << "\n"
387  << " IN_ARG_x_dot_dot = "
388  << inArgs.supports(MEB::IN_ARG_x_dot_dot) << "\n"
389  << " IN_ARG_alpha = "
390  << inArgs.supports(MEB::IN_ARG_alpha) << "\n"
391  << " IN_ARG_beta = "
392  << inArgs.supports(MEB::IN_ARG_beta) << "\n"
393  << " IN_ARG_W_x_dot_dot_coeff = "
394  << inArgs.supports(MEB::IN_ARG_W_x_dot_dot_coeff)
395  << "\n"
396  << " OUT_ARG_f = "
397  << outArgs.supports(MEB::OUT_ARG_f) << "\n"
398  << " OUT_ARG_W = "
399  << outArgs.supports(MEB::OUT_ARG_W) << "\n"
400  << "Implicit Second Order ODE requires:\n"
401  << " IN_ARG_x = true\n"
402  << " IN_ARG_x_dot = true\n"
403  << " IN_ARG_x_dot_dot = true\n"
404  << " IN_ARG_alpha = true\n"
405  << " IN_ARG_beta = true\n"
406  << " IN_ARG_W_x_dot_dot_coeff = true\n"
407  << " OUT_ARG_f = true\n"
408  << " OUT_ARG_W = true\n");
409 
410  return;
411 }
412 
414 {
416  using Teuchos::RCP;
417 
418  // NOX Solver ParameterList
419  RCP<ParameterList> noxPL = Teuchos::parameterList();
420 
421  // Direction ParameterList
422  RCP<ParameterList> directionPL = Teuchos::parameterList();
423  directionPL->set<std::string>("Method", "Newton");
424  RCP<ParameterList> newtonPL = Teuchos::parameterList();
425  newtonPL->set<std::string>("Forcing Term Method", "Constant");
426  newtonPL->set<bool>("Rescue Bad Newton Solve", 1);
427  directionPL->set("Newton", *newtonPL);
428  noxPL->set("Direction", *directionPL);
429 
430  // Line Search ParameterList
431  RCP<ParameterList> lineSearchPL = Teuchos::parameterList();
432  lineSearchPL->set<std::string>("Method", "Full Step");
433  RCP<ParameterList> fullStepPL = Teuchos::parameterList();
434  fullStepPL->set<double>("Full Step", 1);
435  lineSearchPL->set("Full Step", *fullStepPL);
436  noxPL->set("Line Search", *lineSearchPL);
437 
438  noxPL->set<std::string>("Nonlinear Solver", "Line Search Based");
439 
440  // Printing ParameterList
441  RCP<ParameterList> printingPL = Teuchos::parameterList();
442  printingPL->set<int>("Output Precision", 3);
443  printingPL->set<int>("Output Processor", 0);
444  RCP<ParameterList> outputPL = Teuchos::parameterList();
445  outputPL->set<bool>("Error", 1);
446  outputPL->set<bool>("Warning", 1);
447  outputPL->set<bool>("Outer Iteration", 0);
448  outputPL->set<bool>("Parameters", 0);
449  outputPL->set<bool>("Details", 0);
450  outputPL->set<bool>("Linear Solver Details", 1);
451  outputPL->set<bool>("Stepper Iteration", 1);
452  outputPL->set<bool>("Stepper Details", 1);
453  outputPL->set<bool>("Stepper Parameters", 1);
454  printingPL->set("Output Information", *outputPL);
455  noxPL->set("Printing", *printingPL);
456 
457  // Solver Options ParameterList
458  RCP<ParameterList> solverOptionsPL = Teuchos::parameterList();
459  solverOptionsPL->set<std::string>("Status Test Check Type", "Minimal");
460  noxPL->set("Solver Options", *solverOptionsPL);
461 
462  // Status Tests ParameterList
463  RCP<ParameterList> statusTestsPL = Teuchos::parameterList();
464  statusTestsPL->set<std::string>("Test Type", "Combo");
465  statusTestsPL->set<std::string>("Combo Type", "OR");
466  statusTestsPL->set<int>("Number of Tests", 2);
467  RCP<ParameterList> test0PL = Teuchos::parameterList();
468  test0PL->set<std::string>("Test Type", "NormF");
469  test0PL->set<double>("Tolerance", 1e-08);
470  statusTestsPL->set("Test 0", *test0PL);
471  RCP<ParameterList> test1PL = Teuchos::parameterList();
472  test1PL->set<std::string>("Test Type", "MaxIters");
473  test1PL->set<int>("Maximum Iterations", 10);
474  statusTestsPL->set("Test 1", *test1PL);
475  noxPL->set("Status Tests", *statusTestsPL);
476 
477  // Solver ParameterList
478  RCP<ParameterList> solverPL = Teuchos::parameterList("Default Solver");
479  solverPL->set("NOX", *noxPL);
480 
481  return solverPL;
482 }
483 
484 } // namespace Tempus
485 #endif // Tempus_Stepper_impl_hpp
const std::string & name() const
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
bool is_null(const boost::shared_ptr< T > &p)
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.
void validExplicitODE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate that the model supports explicit ODE evaluation, f(x,t) [=xdot].
virtual void initialize()
Initialize after construction and changing input parameters.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasic() const
Add basic parameters to Steppers ParameterList.
void setUseFSALFalseOnly(bool a)
virtual void checkInitialized()
Check initialization, and error out on failure.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void setUseFSALTrueOnly(bool a)
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getStepperX()
Get Stepper x.
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void setStepperValues(const Teuchos::RCP< Teuchos::ParameterList > pl)
Set Stepper member data from ParameterList.
void validSecondOrderExplicitODE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
void validSecondOrderODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getStepperXDotDot()
Get Stepper xDotDot.
void validImplicitODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate ME supports implicit ODE/DAE evaluation, f(xdot,x,t) [= 0].
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Solution state for integrators and steppers.
std::string toString(const T &t)
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getStepperXDot()
Get Stepper xDot.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const