9 #ifndef Tempus_StepperImplicit_impl_hpp
10 #define Tempus_StepperImplicit_impl_hpp
13 #include "NOX_Thyra.H"
17 template <
class Scalar>
26 "Error - Solver is not set!\n");
27 solver_->setModel(wrapperModel_);
29 this->isInitialized_ =
false;
32 template <
class Scalar>
38 int numStates = solutionHistory->getNumStates();
41 numStates < 1, std::logic_error,
42 "Error - setInitialConditions() needs at least one SolutionState\n"
43 " to set the initial condition. Number of States = "
46 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
47 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
48 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
49 if (xDot() == Teuchos::null) xDot = this->getStepperXDot();
51 auto inArgs = this->wrapperModel_->getNominalValues();
53 RCP<Thyra::VectorBase<Scalar> > initialXDot = initialState->getXDot();
57 !((x != Teuchos::null && initialXDot != Teuchos::null) ||
58 (inArgs.get_x() != Teuchos::null &&
59 inArgs.get_x_dot() != Teuchos::null)),
61 "Error - We need to set the initial conditions for x and xDot from\n"
62 " either initialState or appModel_->getNominalValues::InArgs\n"
63 " (but not from a mixture of the two).\n");
68 if (x == Teuchos::null) {
70 (x == Teuchos::null) && (inArgs.get_x() == Teuchos::null),
72 "Error - setInitialConditions needs the ICs from the "
74 " or getNominalValues()!\n");
77 initialState->setX(x);
82 using Teuchos::rcp_const_cast;
83 if (x == Teuchos::null || xDot == Teuchos::null) {
85 (inArgs.get_x() == Teuchos::null) ||
86 (inArgs.get_x_dot() == Teuchos::null),
88 "Error - setInitialConditions() needs the ICs from the initialState\n"
89 " or getNominalValues()!\n");
91 initialState->setX(x);
92 RCP<Thyra::VectorBase<Scalar> > x_dot =
94 initialState->setXDot(x_dot);
99 std::string icConsistency = this->getICConsistency();
100 if (icConsistency ==
"None") {
102 if (initialState->getXDot() == Teuchos::null)
103 Thyra::assign(xDot.ptr(), Scalar(0.0));
106 if (initialState->getXDotDot() == Teuchos::null)
107 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
110 else if (icConsistency ==
"Zero") {
112 Thyra::assign(xDot.ptr(), Scalar(0.0));
114 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
116 else if (icConsistency ==
"App") {
121 x_dot == Teuchos::null, std::logic_error,
122 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
123 " but 'App' returned a null pointer for xDot!\n");
124 Thyra::assign(xDot.ptr(), *x_dot);
128 inArgs.get_x_dot_dot());
130 xDotDot == Teuchos::null, std::logic_error,
131 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
132 " but 'App' returned a null pointer for xDotDot!\n");
133 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
136 else if (icConsistency ==
"Consistent") {
139 const Scalar time = initialState->getTime();
140 const Scalar dt = initialState->getTimeStep();
141 RCP<TimeDerivative<Scalar> > timeDer = Teuchos::null;
142 const Scalar alpha = Scalar(1.0);
143 const Scalar beta = Scalar(0.0);
148 this->solveImplicitODE(x, xDot, time, p);
153 "Error - Solver failed while determining the initial conditions.\n"
159 true, std::logic_error,
160 "Error - setInitialConditions(): 'Consistent' for second-order ODE\n"
161 " has not been implemented.\n");
166 true, std::logic_error,
167 "Error - setInitialConditions() invalid IC consistency, "
168 << icConsistency <<
".\n");
173 initialState->setIsSynced(
true);
176 if (this->getICConsistencyCheck()) {
177 auto f = initialState->getX()->clone_v();
179 const Scalar time = initialState->getTime();
180 const Scalar dt = initialState->getTimeStep();
181 RCP<TimeDerivative<Scalar> > timeDer = Teuchos::null;
182 const Scalar alpha = Scalar(0.0);
183 const Scalar beta = Scalar(0.0);
187 this->evaluateImplicitODE(f, x, xDot, time, p);
189 Scalar normX = Thyra::norm(*x);
190 Scalar reldiff = Scalar(0.0);
191 if (normX == Scalar(0.0))
192 reldiff = Thyra::norm(*f);
194 reldiff = Thyra::norm(*f) / normX;
197 RCP<Teuchos::FancyOStream> out = this->getOStream();
198 Teuchos::OSTab ostab(out, 1,
"StepperImplicit::setInitialConditions()");
200 *out <<
"\n---------------------------------------------------\n"
201 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
202 <<
" Initial condition PASSED consistency check!\n"
203 <<
" (||f(x,xDot,t)||/||x|| = " << reldiff <<
") < "
204 <<
"(eps = " << eps <<
")" << std::endl
205 <<
"---------------------------------------------------\n"
209 *out <<
"\n---------------------------------------------------\n"
210 <<
"Info -- Stepper = " << this->getStepperType() <<
"\n"
211 <<
" Initial condition FAILED consistency check but continuing!\n"
212 <<
" (||f(x,xDot,t)||/||x|| = " << reldiff <<
") > "
213 <<
"(eps = " << eps <<
")" << std::endl
214 <<
" ||f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
215 <<
" ||x|| = " << Thyra::norm(*x) << std::endl
216 <<
"---------------------------------------------------\n"
222 template <
class Scalar>
227 this->setSolverName(
"Default Solver");
228 auto subPL = sublist(solverPL,
"NOX");
229 solver_->setParameterList(subPL);
231 if (wrapperModel_ != Teuchos::null) solver_->setModel(wrapperModel_);
233 this->isInitialized_ =
false;
236 template <
class Scalar>
241 solver == Teuchos::null, std::logic_error,
242 "Error - Can not set the solver to Teuchos::null.\n");
245 if (wrapperModel_ != Teuchos::null) solver_->setModel(wrapperModel_);
247 this->isInitialized_ =
false;
250 template <
class Scalar>
257 wrapperModel_->setForSolve(x, xDot, time, p, y, index);
261 switch (p->evaluationType_) {
263 if (getZeroInitialGuess()) Thyra::assign(x.ptr(), ST::zero());
264 sStatus = (*solver_).solve(&*x);
269 sStatus = (*solver_).solve(&*xDot);
280 template <
class Scalar>
288 MEB::InArgs<Scalar> inArgs = wrapperModel_->createInArgs();
290 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(xDot);
291 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
292 if (inArgs.supports(MEB::IN_ARG_step_size))
293 inArgs.set_step_size(p->timeStepSize_);
294 if (inArgs.supports(MEB::IN_ARG_alpha)) inArgs.set_alpha(Scalar(0.0));
295 if (inArgs.supports(MEB::IN_ARG_beta)) inArgs.set_beta(Scalar(0.0));
297 MEB::OutArgs<Scalar> outArgs = wrapperModel_->createOutArgs();
300 wrapperModel_->setForSolve(x, xDot, time, p);
302 wrapperModel_->evalModel(inArgs, outArgs);
305 template <
class Scalar>
310 out <<
"--- StepperImplicit ---\n";
311 out <<
" wrapperModel_ = " << wrapperModel_ << std::endl;
312 out <<
" solver_ = " << solver_ << std::endl;
313 out <<
" initialGuess_ = " << initialGuess_ << std::endl;
318 template <
class Scalar>
322 bool isValidSetup =
true;
324 if (wrapperModel_->getAppModel() == Teuchos::null) {
325 isValidSetup =
false;
326 out <<
"The application ModelEvaluator is not set!\n";
329 if (wrapperModel_ == Teuchos::null) {
330 isValidSetup =
false;
331 out <<
"The wrapper ModelEvaluator is not set!\n";
334 if (solver_ == Teuchos::null) {
335 isValidSetup =
false;
336 out <<
"The solver is not set!\n";
339 if (solver_ != Teuchos::null) {
340 if (solver_->getModel() == Teuchos::null) {
341 isValidSetup =
false;
342 out <<
"The solver's ModelEvaluator is not set!\n";
349 template <
class Scalar>
353 return this->getValidParametersBasicImplicit();
356 template <
class Scalar>
360 auto pl = this->getValidParametersBasic();
361 pl->template set<std::string>(
"Solver Name", this->getSolverName());
362 pl->template set<bool>(
"Zero Initial Guess", this->getZeroInitialGuess());
363 auto noxSolverPL = this->getSolver()->getParameterList();
364 auto solverPL = Teuchos::parameterList(this->getSolverName());
365 solverPL->set(
"NOX", *noxSolverPL);
366 pl->set(this->getSolverName(), *solverPL);
371 template <
class Scalar>
375 if (pl != Teuchos::null) {
378 this->setStepperValues(pl);
379 this->setZeroInitialGuess(pl->
get<
bool>(
"Zero Initial Guess",
false));
381 this->setStepperSolverValues(pl);
384 template <
class Scalar>
388 if (pl != Teuchos::null) {
390 std::string solverName = pl->
get<std::string>(
"Solver Name");
392 auto solverPL = Teuchos::parameterList();
393 solverPL = Teuchos::sublist(pl, solverName);
394 this->setSolverName(solverName);
396 Teuchos::sublist(solverPL,
"NOX",
true);
397 getSolver()->setParameterList(noxPL);
403 #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)
#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
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=-1)
Solve implicit ODE, f(x, xDot, t, p) = 0.
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
virtual void setDefaultSolver()
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).
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)