Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
16namespace Tempus {
17
18
19template<class Scalar>
21 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
22{
23 validImplicitODE_DAE(appModel);
24 wrapperModel_ =
25 Teuchos::rcp(new WrapperModelEvaluatorBasic<Scalar>(appModel));
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
35template<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.
58 TEUCHOS_TEST_FOR_EXCEPTION(
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)
140 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
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
146 TEUCHOS_TEST_FOR_EXCEPTION(
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);
177 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
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
211template<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
226template<class Scalar>
228 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > solver)
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
240template<class Scalar>
241const Thyra::SolveStatus<Scalar>
243 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & x,
244 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & xDot,
245 const Scalar time,
246 const Teuchos::RCP<ImplicitODEParameters<Scalar> > & p,
247 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & y,
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
266 Thyra::SolveStatus<Scalar> sStatus;
267 typedef Teuchos::ScalarTraits<Scalar> ST;
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 }
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
289template<class Scalar>
290void
292 Teuchos::RCP<Thyra::VectorBase<Scalar> > & f,
293 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & x,
294 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & xDot,
295 const Scalar time,
296 const Teuchos::RCP<ImplicitODEParameters<Scalar> > & p )
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
316template<class Scalar>
317void StepperImplicit<Scalar>::describe(Teuchos::FancyOStream & out,
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
330template<class Scalar>
331bool StepperImplicit<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
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
362template<class Scalar>
363Teuchos::RCP<const Teuchos::ParameterList>
365{
366 return this->getValidParametersBasicImplicit();
367}
368
369
370template<class Scalar>
371Teuchos::RCP<Teuchos::ParameterList>
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
386template<class Scalar>
389 Teuchos::RCP<Teuchos::ParameterList> pl)
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
401template<class Scalar>
403setStepperSolverValues(Teuchos::RCP<Teuchos::ParameterList> pl)
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);
412 Teuchos::RCP<Teuchos::ParameterList> noxPL =
413 Teuchos::sublist(solverPL,"NOX",true);
414 getSolver()->setParameterList(noxPL);
415 }
416 }
417}
418
419
420} // namespace Tempus
421#endif // Tempus_StepperImplicit_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) override
Set the initial conditions and make them consistent.
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicImplicit() const
void setStepperImplicitValues(Teuchos::RCP< Teuchos::ParameterList > pl)
Set StepperImplicit member data from the ParameterList.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
void setStepperSolverValues(Teuchos::RCP< Teuchos::ParameterList > pl)
Set solver from ParameterList.
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver) override
Set solver.
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).
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.
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.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const override
A ModelEvaluator for residual evaluations given a state. This ModelEvaluator takes a state,...
void validImplicitODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate ME supports implicit ODE/DAE evaluation, f(xdot,x,t) [= 0].
@ SECOND_ORDER_ODE
Stepper integrates second-order ODEs.
@ FIRST_ORDER_ODE
Stepper integrates first-order ODEs.
@ SOLVE_FOR_XDOT_CONST_X
Solve for xDot keeping x constant (for ICs).
@ SOLVE_FOR_X
Solve for x and determine xDot from x.
@ EVALUATE_RESIDUAL
Evaluate residual for the implicit ODE.
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.