Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperExplicit_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_StepperExplicit_impl_hpp
10#define Tempus_StepperExplicit_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13
14
15namespace Tempus {
16
17
18template<class Scalar>
20 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
21{
22 validExplicitODE(appModel);
23 appModel_ = appModel;
24
25 inArgs_ = appModel_->getNominalValues();
26 outArgs_ = appModel_->createOutArgs();
27}
28
29
30template<class Scalar>
32 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
33{
34 using Teuchos::RCP;
35
36 int numStates = solutionHistory->getNumStates();
37
38 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
39 "Error - setInitialConditions() needs at least one SolutionState\n"
40 " to set the initial condition. Number of States = " << numStates);
41
42 if (numStates > 1) {
43 RCP<Teuchos::FancyOStream> out = this->getOStream();
44 Teuchos::OSTab ostab(out,1, "StepperExplicit::setInitialConditions()");
45 *out << "Warning -- SolutionHistory has more than one state!\n"
46 << "Setting the initial conditions on the currentState.\n"<<std::endl;
47 }
48
49 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
50 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
51 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
52 if (xDot == Teuchos::null) xDot = this->getStepperXDot();
53
54
55 inArgs_ = appModel_->getNominalValues();
56 if (this->getOrderODE() == SECOND_ORDER_ODE) {
57 RCP<Thyra::VectorBase<Scalar> > initialXDot = initialState->getXDot();
58 // If initialState has x and xDot set, treat them as the initial conditions.
59 // Otherwise use the x and xDot from getNominalValues() as the ICs.
60 TEUCHOS_TEST_FOR_EXCEPTION(
61 !((x != Teuchos::null && initialXDot != Teuchos::null) ||
62 (inArgs_.get_x() != Teuchos::null &&
63 inArgs_.get_x_dot() != Teuchos::null)), std::logic_error,
64 "Error - We need to set the initial conditions for x and xDot from\n"
65 " either initialState or appModel_->getNominalValues::InArgs\n"
66 " (but not from a mixture of the two).\n");
67 }
68
69 if (this->getOrderODE() == FIRST_ORDER_ODE) {
70 if (x == Teuchos::null) {
71 // Use x from inArgs as ICs.
72 TEUCHOS_TEST_FOR_EXCEPTION( (x == Teuchos::null) &&
73 (inArgs_.get_x() == Teuchos::null), std::logic_error,
74 "Error - setInitialConditions needs the ICs from the SolutionHistory\n"
75 " or getNominalValues()!\n");
76
77 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs_.get_x());
78 initialState->setX(x);
79 }
80 }
81 else if (this->getOrderODE() == SECOND_ORDER_ODE) {
82 using Teuchos::rcp_const_cast;
83 // Use the x and x_dot from getNominalValues() as the ICs.
84 if ( initialState->getX() == Teuchos::null ||
85 initialState->getXDot() == Teuchos::null ) {
86 TEUCHOS_TEST_FOR_EXCEPTION( (inArgs_.get_x() == Teuchos::null) ||
87 (inArgs_.get_x_dot() == Teuchos::null), std::logic_error,
88 "Error - setInitialConditions() needs the ICs from the initialState\n"
89 " or getNominalValues()!\n");
90 x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs_.get_x());
91 initialState->setX(x);
92 RCP<Thyra::VectorBase<Scalar> > x_dot
93 = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs_.get_x_dot());
94 initialState->setXDot(x_dot);
95 }
96 }
97
98 // Perform IC Consistency
99 std::string icConsistency = this->getICConsistency();
100 if (icConsistency == "None") {
101 if (this->getOrderODE() == FIRST_ORDER_ODE) {
102 if (initialState->getXDot() == Teuchos::null) {
103 RCP<Teuchos::FancyOStream> out = this->getOStream();
104 Teuchos::OSTab ostab(out,1,"StepperExplicit::setInitialConditions()");
105 *out << "Warning -- Requested IC consistency of 'None' but\n"
106 << " initialState does not have an xDot.\n"
107 << " Setting a 'Consistent' xDot!\n" << std::endl;
108 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
109 evaluateExplicitODE(xDot, x, initialState->getTime(), p);
110 initialState->setIsSynced(true);
111 }
112 }
113 else if (this->getOrderODE() == SECOND_ORDER_ODE) {
114 if (initialState->getXDotDot() == Teuchos::null) {
115 RCP<Teuchos::FancyOStream> out = this->getOStream();
116 Teuchos::OSTab ostab(out,1,"StepperExplicit::setInitialConditions()");
117 *out << "Warning -- Requested IC consistency of 'None' but\n"
118 << " initialState does not have an xDotDot.\n"
119 << " Setting a 'Consistent' xDotDot!\n" << std::endl;
120 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
121 this->evaluateExplicitODE(this->getStepperXDotDot(initialState), x,
122 initialState->getXDot(),
123 initialState->getTime(), p);
124 initialState->setIsSynced(true);
125 }
126 }
127 }
128 else if (icConsistency == "Zero") {
129 if (this->getOrderODE() == FIRST_ORDER_ODE)
130 Thyra::assign(xDot.ptr(), Scalar(0.0));
131 else if (this->getOrderODE() == SECOND_ORDER_ODE)
132 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
133 }
134 else if (icConsistency == "App") {
135 if (this->getOrderODE() == FIRST_ORDER_ODE) {
136 auto x_dot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
137 inArgs_.get_x_dot());
138 TEUCHOS_TEST_FOR_EXCEPTION(xDot == Teuchos::null, std::logic_error,
139 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
140 " but 'App' returned a null pointer for xDot!\n");
141 Thyra::assign(xDot.ptr(), *x_dot);
142 }
143 else if (this->getOrderODE() == SECOND_ORDER_ODE) {
144 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
145 inArgs_.get_x_dot_dot());
146 TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
147 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
148 " but 'App' returned a null pointer for xDotDot!\n");
149 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
150 }
151 }
152 else if (icConsistency == "Consistent") {
153 if (this->getOrderODE() == FIRST_ORDER_ODE) {
154 // Evaluate xDot = f(x,t).
155 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
156 evaluateExplicitODE(xDot, x, initialState->getTime(), p);
157 }
158 else if (this->getOrderODE() == SECOND_ORDER_ODE) {
159 // Evaluate xDotDot = f(x,xDot,t).
160 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
161 this->evaluateExplicitODE(initialState->getXDotDot(), x,
162 initialState->getXDot(),
163 initialState->getTime(), p);
164 }
165 }
166 else {
167 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
168 "Error - setInitialConditions() invalid IC consistency, '"
169 << icConsistency << "'.\n");
170 }
171
172 // At this point, x, and xDot (and xDotDot) sync'ed or consistent
173 // at the same time level for the initialState.
174 initialState->setIsSynced(true);
175
176 // Test for consistency.
177 if (this->getICConsistencyCheck()) {
178 if (this->getOrderODE() == FIRST_ORDER_ODE) {
179 auto f = initialState->getX()->clone_v();
180 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
181 evaluateExplicitODE(f, x, initialState->getTime(), p);
182 Thyra::Vp_StV(f.ptr(), Scalar(-1.0), *(xDot));
183 Scalar normX = Thyra::norm(*x);
184 Scalar reldiff = Scalar(0.0);
185 if (normX == Scalar(0.0)) reldiff = Thyra::norm(*f);
186 else reldiff = Thyra::norm(*f)/normX;
187
188 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
189 RCP<Teuchos::FancyOStream> out = this->getOStream();
190 Teuchos::OSTab ostab(out,1,"StepperExplicit::setInitialConditions()");
191 if (reldiff < eps) {
192 *out << "\n---------------------------------------------------\n"
193 << "Info -- Stepper = " << this->getStepperType() << "\n"
194 << " Initial condition PASSED consistency check!\n"
195 << " (||xDot-f(x,t)||/||x|| = " << reldiff << ") < "
196 << "(eps = " << eps << ")" << std::endl
197 << "---------------------------------------------------\n"<<std::endl;
198 } else {
199 *out << "\n---------------------------------------------------\n"
200 << "Info -- Stepper = " << this->getStepperType() << "\n"
201 << " Initial condition FAILED consistency check but continuing!\n"
202 << " (||xDot-f(x,t)||/||x|| = " << reldiff << ") > "
203 << "(eps = " << eps << ")" << std::endl
204 << " ||xDot-f(x,t)|| = " << Thyra::norm(*f) << std::endl
205 << " ||x|| = " << Thyra::norm(*x) << std::endl
206 << "---------------------------------------------------\n"<<std::endl;
207 }
208 }
209 else if (this->getOrderODE() == SECOND_ORDER_ODE) {
210 auto xDotDot = initialState->getXDotDot();
211 auto f = initialState->getX()->clone_v();
212 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(0.0));
213 this->evaluateExplicitODE(f, x, initialState->getXDot(),
214 initialState->getTime(), p);
215 Thyra::Vp_StV(f.ptr(), Scalar(-1.0), *(xDotDot));
216 Scalar normX = Thyra::norm(*x);
217 Scalar reldiff = Scalar(0.0);
218 if (normX == Scalar(0.0)) reldiff = Thyra::norm(*f);
219 else reldiff = Thyra::norm(*f)/normX;
220
221 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
222 RCP<Teuchos::FancyOStream> out = this->getOStream();
223 Teuchos::OSTab ostab(out,1,"StepperExplicit::setInitialConditions()");
224 if (reldiff < eps) {
225 *out << "\n---------------------------------------------------\n"
226 << "Info -- Stepper = " << this->getStepperType() << "\n"
227 << " Initial condition PASSED consistency check!\n"
228 << " (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff << ") > "
229 << "(eps = " << eps << ")" << std::endl
230 << "---------------------------------------------------\n"<<std::endl;
231 } else {
232 *out << "\n---------------------------------------------------\n"
233 << "Info -- Stepper = " << this->getStepperType() << "\n"
234 << "Initial condition FAILED consistency check but continuing!\n"
235 << " (||xDotDot-f(x,xDot,t)||/||x|| = " << reldiff << ") > "
236 << "(eps = " << eps << ")" << std::endl
237 << " ||xDotDot-f(x,xDot,t)|| = " << Thyra::norm(*f) << std::endl
238 << " ||x|| = " << Thyra::norm(*x) << std::endl
239 << "---------------------------------------------------\n"<<std::endl;
240 }
241 }
242 }
243}
244
245template<class Scalar>
247 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > /* solver */)
248{
249 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
250 Teuchos::OSTab ostab(out,1,"StepperExplicit::setSolver()");
251 *out << "Warning -- No solver to set for StepperExplicit "
252 << "(i.e., explicit method).\n" << std::endl;
253 return;
254}
255template<class Scalar>
256void
259 Teuchos::RCP<const Thyra::VectorBase<Scalar> > x,
260 const Scalar time,
261 const Teuchos::RCP<ExplicitODEParameters<Scalar> > & p )
262{
263 typedef Thyra::ModelEvaluatorBase MEB;
264
265 inArgs_.set_x(x);
266 if (inArgs_.supports(MEB::IN_ARG_t)) inArgs_.set_t(time);
267 if (inArgs_.supports(MEB::IN_ARG_step_size))
268 inArgs_.set_step_size(p->timeStepSize_);
269 if (inArgs_.supports(MEB::IN_ARG_stage_number))
270 inArgs_.set_stage_number(p->stageNumber_);
271
272 // For model evaluators whose state function f(x, xDot, t) describes
273 // an implicit ODE, and which accept an optional xDot input argument,
274 // make sure the latter is set to null in order to request the evaluation
275 // of a state function corresponding to the explicit ODE formulation
276 // xDot = f(x, t).
277 if (inArgs_.supports(MEB::IN_ARG_x_dot)) inArgs_.set_x_dot(Teuchos::null);
278
279 outArgs_.set_f(xDot);
280
281 appModel_->evalModel(inArgs_, outArgs_);
282}
283
284template<class Scalar>
285void
287evaluateExplicitODE(Teuchos::RCP< Thyra::VectorBase<Scalar> > xDotDot,
288 Teuchos::RCP<const Thyra::VectorBase<Scalar> > x,
289 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDot,
290 const Scalar time,
291 const Teuchos::RCP<ExplicitODEParameters<Scalar> > & p )
292{
293 typedef Thyra::ModelEvaluatorBase MEB;
294
295 inArgs_.set_x(x);
296 if (inArgs_.supports(MEB::IN_ARG_x_dot)) inArgs_.set_x_dot(xDot);
297 if (inArgs_.supports(MEB::IN_ARG_t)) inArgs_.set_t(time);
298 if (inArgs_.supports(MEB::IN_ARG_step_size))
299 inArgs_.set_step_size(p->timeStepSize_);
300 if (inArgs_.supports(MEB::IN_ARG_stage_number))
301 inArgs_.set_stage_number(p->stageNumber_);
302
303 // For model evaluators whose state function f(x, xDot, xDotDot, t) describes
304 // an implicit ODE, and which accept an optional xDotDot input argument,
305 // make sure the latter is set to null in order to request the evaluation
306 // of a state function corresponding to the explicit ODE formulation
307 // xDotDot = f(x, xDot, t).
308 if (inArgs_.supports(MEB::IN_ARG_x_dot_dot))
309 inArgs_.set_x_dot_dot(Teuchos::null);
310
311 outArgs_.set_f(xDotDot);
312
313 appModel_->evalModel(inArgs_, outArgs_);
314}
315
316
317
318template<class Scalar>
319void StepperExplicit<Scalar>::describe(Teuchos::FancyOStream & out,
320 const Teuchos::EVerbosityLevel verbLevel) const
321{
322 auto l_out = Teuchos::fancyOStream( out.getOStream() );
323 Teuchos::OSTab ostab(*l_out, 2, this->description());
324 l_out->setOutputToRootOnly(0);
325
326 *l_out << "--- StepperExplicit ---\n"
327 << " appModel_ = " << appModel_ << std::endl
328 << " inArgs_ = " << inArgs_ << std::endl
329 << " outArgs_ = " << outArgs_ << std::endl;
330}
331
332
333template<class Scalar>
334bool StepperExplicit<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
335{
336 out.setOutputToRootOnly(0);
337 bool isValidSetup = true;
338
339 if (appModel_ == Teuchos::null) {
340 isValidSetup = false;
341 out << "The application ModelEvaluator is not set!\n";
342 }
343
344 return isValidSetup;
345}
346
347
348template<class Scalar>
350setStepperExplicitValues(Teuchos::RCP<Teuchos::ParameterList> pl)
351{
352 if (pl != Teuchos::null) {
353 pl->validateParametersAndSetDefaults(*this->getValidParameters(), 0);
354 this->setStepperValues(pl);
355 }
356}
357
358
359} // namespace Tempus
360#endif // Tempus_StepperExplicit_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
virtual void evaluateExplicitODE(Teuchos::RCP< Thyra::VectorBase< Scalar > > xDot, Teuchos::RCP< const Thyra::VectorBase< Scalar > > x, const Scalar time, const Teuchos::RCP< ExplicitODEParameters< Scalar > > &p)
Evaluate xDot = f(x,t).
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions, make them consistent, and set needed memory.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set model.
void setStepperExplicitValues(Teuchos::RCP< Teuchos::ParameterList > pl)
Set StepperExplicit member data from the ParameterList.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void setSolver(Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > solver)
Set solver.
@ SECOND_ORDER_ODE
Stepper integrates second-order ODEs.
@ FIRST_ORDER_ODE
Stepper integrates first-order ODEs.
void validExplicitODE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate that the model supports explicit ODE evaluation, f(x,t) [=xdot].