9 #ifndef Tempus_StepperImplicit_impl_hpp
10 #define Tempus_StepperImplicit_impl_hpp
13 #include "NOX_Thyra.H"
19 template<
class Scalar>
28 "Error - Solver is not set!\n");
29 solver_->setModel(wrapperModel_);
31 this->isInitialized_ =
false;
35 template<
class Scalar>
41 int numStates = solutionHistory->getNumStates();
44 "Error - setInitialConditions() needs at least one SolutionState\n"
45 " to set the initial condition. Number of States = " << numStates);
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();
53 auto inArgs = this->wrapperModel_->getNominalValues();
55 RCP<Thyra::VectorBase<Scalar> > initialXDot = initialState->getXDot();
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");
69 if (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");
76 initialState->setX(x);
81 using Teuchos::rcp_const_cast;
82 if ( x == Teuchos::null || xDot == 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");
88 initialState->setX(x);
89 RCP<Thyra::VectorBase<Scalar> > x_dot =
91 initialState->setXDot(x_dot);
97 std::string icConsistency = this->getICConsistency();
98 if (icConsistency ==
"None") {
100 if (initialState->getXDot() == Teuchos::null)
101 Thyra::assign(xDot.ptr(), Scalar(0.0));
104 if (initialState->getXDotDot() == Teuchos::null)
105 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
108 else if (icConsistency ==
"Zero") {
110 Thyra::assign(xDot.ptr(), Scalar(0.0));
112 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
114 else if (icConsistency ==
"App") {
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);
125 inArgs.get_x_dot_dot());
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);
132 else if (icConsistency ==
"Consistent") {
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);
139 const Scalar beta = Scalar(0.0);
144 this->solveImplicitODE(x, xDot, time, p);
148 "Error - Solver failed while determining the initial conditions.\n"
154 "Error - setInitialConditions(): 'Consistent' for second-order ODE\n"
155 " has not been implemented.\n");
160 "Error - setInitialConditions() invalid IC consistency, "
161 << icConsistency <<
".\n");
166 initialState->setIsSynced(
true);
169 if (this->getICConsistencyCheck()) {
170 auto f = initialState->getX()->clone_v();
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);
180 this->evaluateImplicitODE(f, x, xDot, time, p);
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;
188 RCP<Teuchos::FancyOStream> out = this->getOStream();
189 Teuchos::OSTab ostab(out,1,
"StepperImplicit::setInitialConditions()");
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;
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;
211 template<
class Scalar>
216 this->setSolverName(
"Default Solver");
217 auto subPL = sublist(solverPL,
"NOX");
218 solver_->setParameterList(subPL);
220 if (wrapperModel_ != Teuchos::null) solver_->setModel(wrapperModel_);
222 this->isInitialized_ =
false;
226 template<
class Scalar>
231 "Error - Can not set the solver to Teuchos::null.\n");
234 if (wrapperModel_ != Teuchos::null) solver_->setModel(wrapperModel_);
236 this->isInitialized_ =
false;
240 template<
class Scalar>
251 MEB::InArgs<Scalar> inArgs = wrapperModel_->getInArgs();
252 MEB::OutArgs<Scalar> outArgs = wrapperModel_->getOutArgs();
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_);
264 wrapperModel_->setForSolve(p->timeDer_, inArgs, outArgs, p->evaluationType_);
268 switch (p->evaluationType_)
271 if (getZeroInitialGuess()) Thyra::assign(x.ptr(), ST::zero());
272 sStatus = (*solver_).solve(&*x);
277 sStatus = (*solver_).solve(&*xDot);
289 template<
class Scalar>
299 MEB::InArgs<Scalar> inArgs = wrapperModel_->getInArgs();
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));
307 MEB::OutArgs<Scalar> outArgs = wrapperModel_->getOutArgs();
310 wrapperModel_->setForSolve(Teuchos::null,inArgs,outArgs,p->evaluationType_);
312 wrapperModel_->evalModel(inArgs, outArgs);
316 template<
class Scalar>
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_ = "
330 template<
class Scalar>
334 bool isValidSetup =
true;
336 if (wrapperModel_->getAppModel() == Teuchos::null) {
337 isValidSetup =
false;
338 out <<
"The application ModelEvaluator is not set!\n";
341 if (wrapperModel_ == Teuchos::null) {
342 isValidSetup =
false;
343 out <<
"The wrapper ModelEvaluator is not set!\n";
346 if (solver_ == Teuchos::null) {
347 isValidSetup =
false;
348 out <<
"The solver is not set!\n";
351 if (solver_ != Teuchos::null) {
352 if (solver_->getModel() == Teuchos::null) {
353 isValidSetup =
false;
354 out <<
"The solver's ModelEvaluator is not set!\n";
362 template<
class Scalar>
366 return this->getValidParametersBasicImplicit();
370 template<
class Scalar>
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);
386 template<
class Scalar>
391 if (pl != Teuchos::null) {
394 this->setStepperValues(pl);
395 this->setZeroInitialGuess(pl->
get<
bool>(
"Zero Initial Guess",
false));
397 this->setStepperSolverValues(pl);
401 template<
class Scalar>
405 if (pl != Teuchos::null) {
407 std::string solverName = pl->
get<std::string>(
"Solver Name");
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);
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
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)