Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_StepperImplicit_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_StepperImplicit_impl_hpp
10 #define Tempus_StepperImplicit_impl_hpp
11 
12 // Thrya
13 #include "NOX_Thyra.H"
14 
15 
16 namespace Tempus {
17 
18 
19 template<class Scalar>
21  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
22 {
23  validImplicitODE_DAE(appModel);
24  wrapperModel_ =
26 
27  TEUCHOS_TEST_FOR_EXCEPTION(solver_ == Teuchos::null, std::logic_error,
28  "Error - Solver is not set!\n");
29  solver_->setModel(wrapperModel_);
30 
31  this->isInitialized_ = false;
32 }
33 
34 
35 template<class Scalar>
37  const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
38 {
39  using Teuchos::RCP;
40 
41  int numStates = solutionHistory->getNumStates();
42 
43  TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
44  "Error - setInitialConditions() needs at least one SolutionState\n"
45  " to set the initial condition. Number of States = " << numStates);
46 
47  RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
48  RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
49  RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
50  if (xDot() == Teuchos::null) xDot = this->getStepperXDot();
51 
52 
53  auto inArgs = this->wrapperModel_->getNominalValues();
54  if (this->getOrderODE() == SECOND_ORDER_ODE) {
55  RCP<Thyra::VectorBase<Scalar> > initialXDot = initialState->getXDot();
56  // If initialState has x and xDot set, treat them as the initial conditions.
57  // Otherwise use the x and xDot from getNominalValues() as the ICs.
59  !((x != Teuchos::null && initialXDot != Teuchos::null) ||
60  (inArgs.get_x() != Teuchos::null &&
61  inArgs.get_x_dot() != Teuchos::null)), std::logic_error,
62  "Error - We need to set the initial conditions for x and xDot from\n"
63  " either initialState or appModel_->getNominalValues::InArgs\n"
64  " (but not from a mixture of the two).\n");
65  }
66 
67  if (this->getOrderODE() == FIRST_ORDER_ODE) {
68  // Use x from inArgs as ICs.
69  if (x == Teuchos::null) {
70  TEUCHOS_TEST_FOR_EXCEPTION( (x == Teuchos::null) &&
71  (inArgs.get_x() == Teuchos::null), std::logic_error,
72  "Error - setInitialConditions needs the ICs from the SolutionHistory\n"
73  " or getNominalValues()!\n");
74 
75  x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
76  initialState->setX(x);
77  }
78  }
79  else if (this->getOrderODE() == SECOND_ORDER_ODE) {
80  // Use the x and xDot from getNominalValues() as the ICs.
81  using Teuchos::rcp_const_cast;
82  if ( x == Teuchos::null || xDot == Teuchos::null ) {
83  TEUCHOS_TEST_FOR_EXCEPTION( (inArgs.get_x() == Teuchos::null) ||
84  (inArgs.get_x_dot() == Teuchos::null), std::logic_error,
85  "Error - setInitialConditions() needs the ICs from the initialState\n"
86  " or getNominalValues()!\n");
87  x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
88  initialState->setX(x);
89  RCP<Thyra::VectorBase<Scalar> > x_dot =
90  rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
91  initialState->setXDot(x_dot);
92  }
93  }
94 
95 
96  // Perform IC Consistency
97  std::string icConsistency = this->getICConsistency();
98  if (icConsistency == "None") {
99  if (this->getOrderODE() == FIRST_ORDER_ODE) {
100  if (initialState->getXDot() == Teuchos::null)
101  Thyra::assign(xDot.ptr(), Scalar(0.0));
102  }
103  else if (this->getOrderODE() == SECOND_ORDER_ODE) {
104  if (initialState->getXDotDot() == Teuchos::null)
105  Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
106  }
107  }
108  else if (icConsistency == "Zero") {
109  if (this->getOrderODE() == FIRST_ORDER_ODE)
110  Thyra::assign(xDot.ptr(), Scalar(0.0));
111  else if (this->getOrderODE() == SECOND_ORDER_ODE)
112  Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
113  }
114  else if (icConsistency == "App") {
115  if (this->getOrderODE() == FIRST_ORDER_ODE) {
116  auto x_dot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
117  inArgs.get_x_dot());
118  TEUCHOS_TEST_FOR_EXCEPTION(x_dot == Teuchos::null, std::logic_error,
119  "Error - setInitialConditions() requested 'App' for IC consistency,\n"
120  " but 'App' returned a null pointer for xDot!\n");
121  Thyra::assign(xDot.ptr(), *x_dot);
122  }
123  else if (this->getOrderODE() == SECOND_ORDER_ODE) {
124  auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
125  inArgs.get_x_dot_dot());
126  TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
127  "Error - setInitialConditions() requested 'App' for IC consistency,\n"
128  " but 'App' returned a null pointer for xDotDot!\n");
129  Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
130  }
131  }
132  else if (icConsistency == "Consistent") {
133  if (this->getOrderODE() == FIRST_ORDER_ODE) {
134  // Solve f(x, xDot,t) = 0.
135  const Scalar time = initialState->getTime();
136  const Scalar dt = initialState->getTimeStep();
137  RCP<TimeDerivative<Scalar> > timeDer = Teuchos::null;
138  const Scalar alpha = Scalar(1.0); // d(xDot)/d(xDot)
139  const Scalar beta = Scalar(0.0); // d(x )/d(xDot)
141  timeDer, dt, alpha, beta, SOLVE_FOR_XDOT_CONST_X));
142 
143  const Thyra::SolveStatus<Scalar> sStatus =
144  this->solveImplicitODE(x, xDot, time, p);
145 
147  sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED, std::logic_error,
148  "Error - Solver failed while determining the initial conditions.\n"
149  " Solver status is "<<Thyra::toString(sStatus.solveStatus)<<".\n");
150  }
151  else if (this->getOrderODE() == SECOND_ORDER_ODE) {
152 
153  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
154  "Error - setInitialConditions(): 'Consistent' for second-order ODE\n"
155  " has not been implemented.\n");
156  }
157  }
158  else {
159  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
160  "Error - setInitialConditions() invalid IC consistency, "
161  << icConsistency << ".\n");
162  }
163 
164  // At this point, x and xDot are sync'ed or consistent
165  // at the same time level for the initialState.
166  initialState->setIsSynced(true);
167 
168  // Test for consistency.
169  if (this->getICConsistencyCheck()) {
170  auto f = initialState->getX()->clone_v();
171 
172  const Scalar time = initialState->getTime();
173  const Scalar dt = initialState->getTimeStep();
174  RCP<TimeDerivative<Scalar> > timeDer = Teuchos::null;
175  const Scalar alpha = Scalar(0.0);
176  const Scalar beta = Scalar(0.0);
178  timeDer, dt, alpha, beta, EVALUATE_RESIDUAL));
179 
180  this->evaluateImplicitODE(f, x, xDot, time, p);
181 
182  Scalar normX = Thyra::norm(*x);
183  Scalar reldiff = Scalar(0.0);
184  if (normX == Scalar(0.0)) reldiff = Thyra::norm(*f);
185  else reldiff = Thyra::norm(*f)/normX;
186 
187  Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
188  RCP<Teuchos::FancyOStream> out = this->getOStream();
189  Teuchos::OSTab ostab(out,1,"StepperImplicit::setInitialConditions()");
190  if (reldiff < eps) {
191  *out << "\n---------------------------------------------------\n"
192  << "Info -- Stepper = " << this->getStepperType() << "\n"
193  << " Initial condition PASSED consistency check!\n"
194  << " (||f(x,xDot,t)||/||x|| = " << reldiff << ") < "
195  << "(eps = " << eps << ")" << std::endl
196  << "---------------------------------------------------\n"<<std::endl;
197  } else {
198  *out << "\n---------------------------------------------------\n"
199  << "Info -- Stepper = " << this->getStepperType() << "\n"
200  << " Initial condition FAILED consistency check but continuing!\n"
201  << " (||f(x,xDot,t)||/||x|| = " << reldiff << ") > "
202  << "(eps = " << eps << ")" << std::endl
203  << " ||f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
204  << " ||x|| = " << Thyra::norm(*x) << std::endl
205  << "---------------------------------------------------\n"<<std::endl;
206  }
207  }
208 }
209 
210 
211 template<class Scalar>
213 {
214  solver_ = rcp(new Thyra::NOXNonlinearSolver());
215  auto solverPL = defaultSolverParameters();
216  this->setSolverName("Default Solver");
217  auto subPL = sublist(solverPL, "NOX");
218  solver_->setParameterList(subPL);
219 
220  if (wrapperModel_ != Teuchos::null) solver_->setModel(wrapperModel_);
221 
222  this->isInitialized_ = false;
223 }
224 
225 
226 template<class Scalar>
229 {
230  TEUCHOS_TEST_FOR_EXCEPTION(solver == Teuchos::null, std::logic_error,
231  "Error - Can not set the solver to Teuchos::null.\n");
232 
233  solver_ = solver;
234  if (wrapperModel_ != Teuchos::null) solver_->setModel(wrapperModel_);
235 
236  this->isInitialized_ = false;
237 }
238 
239 
240 template<class Scalar>
245  const Scalar time,
248  const int index )
249 {
250  typedef Thyra::ModelEvaluatorBase MEB;
251  MEB::InArgs<Scalar> inArgs = wrapperModel_->getInArgs();
252  MEB::OutArgs<Scalar> outArgs = wrapperModel_->getOutArgs();
253  inArgs.set_x(x);
254  if ( y != Teuchos::null ) inArgs.set_p(index, y);
255  if (inArgs.supports(MEB::IN_ARG_x_dot )) inArgs.set_x_dot (xDot);
256  if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (time);
257  if (inArgs.supports(MEB::IN_ARG_step_size))
258  inArgs.set_step_size(p->timeStepSize_);
259  if (inArgs.supports(MEB::IN_ARG_alpha )) inArgs.set_alpha (p->alpha_);
260  if (inArgs.supports(MEB::IN_ARG_beta )) inArgs.set_beta (p->beta_);
261  if (inArgs.supports(MEB::IN_ARG_stage_number))
262  inArgs.set_stage_number(p->stageNumber_);
263 
264  wrapperModel_->setForSolve(p->timeDer_, inArgs, outArgs, p->evaluationType_);
265 
268  switch (p->evaluationType_)
269  {
270  case SOLVE_FOR_X: {
271  if (getZeroInitialGuess()) Thyra::assign(x.ptr(), ST::zero());
272  sStatus = (*solver_).solve(&*x);
273  break;
274  }
275  case SOLVE_FOR_XDOT_CONST_X: {
276  //if (getZeroInitialGuess()) Thyra::assign(xDot.ptr(), ST::zero());
277  sStatus = (*solver_).solve(&*xDot);
278  break;
279  }
280  default: {
281  TEUCHOS_TEST_FOR_EXCEPT("Invalid EVALUATION_TYPE!");
282  }
283  }
284 
285  return sStatus;
286 }
287 
288 
289 template<class Scalar>
290 void
295  const Scalar time,
297 {
298  typedef Thyra::ModelEvaluatorBase MEB;
299  MEB::InArgs<Scalar> inArgs = wrapperModel_->getInArgs();
300  inArgs.set_x(x);
301  if (inArgs.supports(MEB::IN_ARG_x_dot )) inArgs.set_x_dot (xDot);
302  if (inArgs.supports(MEB::IN_ARG_t )) inArgs.set_t (time);
303  if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(p->timeStepSize_);
304  if (inArgs.supports(MEB::IN_ARG_alpha )) inArgs.set_alpha (Scalar(0.0));
305  if (inArgs.supports(MEB::IN_ARG_beta )) inArgs.set_beta (Scalar(0.0));
306 
307  MEB::OutArgs<Scalar> outArgs = wrapperModel_->getOutArgs();
308  outArgs.set_f(f);
309 
310  wrapperModel_->setForSolve(Teuchos::null,inArgs,outArgs,p->evaluationType_);
311 
312  wrapperModel_->evalModel(inArgs, outArgs);
313 }
314 
315 
316 template<class Scalar>
318  const Teuchos::EVerbosityLevel verbLevel) const
319 {
320  out.setOutputToRootOnly(0);
321  out << "--- StepperImplicit ---\n";
322  out << " wrapperModel_ = " << wrapperModel_ << std::endl;
323  out << " solver_ = " << solver_ << std::endl;
324  out << " initialGuess_ = " << initialGuess_ << std::endl;
325  out << " zeroInitialGuess_ = "
326  << Teuchos::toString(zeroInitialGuess_) << std::endl;
327 }
328 
329 
330 template<class Scalar>
332 {
333  out.setOutputToRootOnly(0);
334  bool isValidSetup = true;
335 
336  if (wrapperModel_->getAppModel() == Teuchos::null) {
337  isValidSetup = false;
338  out << "The application ModelEvaluator is not set!\n";
339  }
340 
341  if (wrapperModel_ == Teuchos::null) {
342  isValidSetup = false;
343  out << "The wrapper ModelEvaluator is not set!\n";
344  }
345 
346  if (solver_ == Teuchos::null) {
347  isValidSetup = false;
348  out << "The solver is not set!\n";
349  }
350 
351  if (solver_ != Teuchos::null) {
352  if (solver_->getModel() == Teuchos::null) {
353  isValidSetup = false;
354  out << "The solver's ModelEvaluator is not set!\n";
355  }
356  }
357 
358  return isValidSetup;
359 }
360 
361 
362 template<class Scalar>
365 {
366  return this->getValidParametersBasicImplicit();
367 }
368 
369 
370 template<class Scalar>
373 {
374  auto pl = this->getValidParametersBasic();
375  pl->template set<std::string>("Solver Name", this->getSolverName());
376  pl->template set<bool>("Zero Initial Guess", this->getZeroInitialGuess());
377  auto noxSolverPL = this->getSolver()->getParameterList();
378  auto solverPL = Teuchos::parameterList(this->getSolverName());
379  solverPL->set("NOX", *noxSolverPL);
380  pl->set(this->getSolverName(), *solverPL);
381 
382  return pl;
383 }
384 
385 
386 template<class Scalar>
390 {
391  if (pl != Teuchos::null) {
392  // Can not validate because of optional Parameters, e.g., 'Solver Name'.
393  //pl->validateParametersAndSetDefaults(*this->getValidParameters());
394  this->setStepperValues(pl);
395  this->setZeroInitialGuess(pl->get<bool>("Zero Initial Guess", false));
396  }
397  this->setStepperSolverValues(pl);
398 }
399 
400 
401 template<class Scalar>
404 {
405  if (pl != Teuchos::null) {
406  setDefaultSolver();
407  std::string solverName = pl->get<std::string>("Solver Name");
408  if ( pl->isSublist(solverName) ) {
409  auto solverPL = Teuchos::parameterList();
410  solverPL = Teuchos::sublist(pl, solverName);
411  this->setSolverName(solverName);
413  Teuchos::sublist(solverPL,"NOX",true);
414  getSolver()->setParameterList(noxPL);
415  }
416  }
417 }
418 
419 
420 } // namespace Tempus
421 #endif // Tempus_StepperImplicit_impl_hpp
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver) override
Set solver.
Evaluate residual for the implicit ODE.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const override
T & get(const std::string &name, T def_value)
const Thyra::SolveStatus< Scalar > solveImplicitODE(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xDot, const Scalar time, const Teuchos::RCP< ImplicitODEParameters< Scalar > > &p, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &y=Teuchos::null, const int index=0)
Solve implicit ODE, f(x, xDot, t, p) = 0.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setStepperImplicitValues(Teuchos::RCP< Teuchos::ParameterList > pl)
Set StepperImplicit member data from the ParameterList.
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel) override
Set the model.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
void evaluateImplicitODE(Teuchos::RCP< Thyra::VectorBase< Scalar > > &f, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xDot, const Scalar time, const Teuchos::RCP< ImplicitODEParameters< Scalar > > &p)
Evaluate implicit ODE residual, f(x, xDot, t, p).
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Stepper integrates second-order ODEs.
bool isSublist(const std::string &name) const
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
void setStepperSolverValues(Teuchos::RCP< Teuchos::ParameterList > pl)
Set solver from ParameterList.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Solve for xDot keeping x constant (for ICs).
ESolveStatus solveStatus
Stepper integrates first-order ODEs.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicImplicit() const
void validImplicitODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate ME supports implicit ODE/DAE evaluation, f(xdot,x,t) [= 0].
A ModelEvaluator for residual evaluations given a state. This ModelEvaluator takes a state...
Solve for x and determine xDot from x.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string toString(const T &t)