Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperExplicitRK_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_StepperExplicitRK_impl_hpp
10#define Tempus_StepperExplicitRK_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13
14
16
17
18namespace Tempus {
19
20
21template<class Scalar>
23{
24 this->setUseEmbedded(false);
25 this->setStageNumber(-1);
26 this->setAppAction(Teuchos::null);
27}
28
29
30template<class Scalar>
32 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
33 bool useFSAL,
34 std::string ICConsistency,
35 bool ICConsistencyCheck,
36 bool useEmbedded,
37 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction)
38{
39 this->setUseFSAL( useFSAL);
40 this->setICConsistency( ICConsistency);
41 this->setICConsistencyCheck( ICConsistencyCheck);
42 this->setUseEmbedded( useEmbedded);
43 this->setStageNumber(-1);
44 this->setErrorNorm();
45
46 this->setAppAction(stepperRKAppAction);
47
48 if (appModel != Teuchos::null) {
49 this->setModel(appModel);
50 this->initialize();
51 }
52}
53
54
55template<class Scalar>
57 const Teuchos::RCP<SolutionHistory<Scalar> >& sh) const
58{
59
60 Scalar dt = Scalar(1.0e+99);
61 if (!this->getUseEmbedded()) return dt;
62
63 Teuchos::RCP<SolutionState<Scalar> > currentState=sh->getCurrentState();
64 const int order = currentState->getOrder();
65 const Scalar time = currentState->getTime();
66 const Scalar errorRel = currentState->getTolRel();
67 const Scalar errorAbs = currentState->getTolAbs();
68
69 Teuchos::RCP<Thyra::VectorBase<Scalar> > stageX, scratchX;
70 stageX = Thyra::createMember(this->appModel_->get_f_space());
71 scratchX = Thyra::createMember(this->appModel_->get_f_space());
72 Thyra::assign(stageX.ptr(), *(currentState->getX()));
73
74 std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar> > > stageXDot(2);
75 for (int i=0; i<2; ++i) {
76 stageXDot[i] = Thyra::createMember(this->appModel_->get_f_space());
77 assign(stageXDot[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
78 }
79
80 // A: one function evaluation at F(t_0, X_0)
81 typedef Thyra::ModelEvaluatorBase MEB;
82 MEB::InArgs<Scalar> inArgs = this->appModel_->getNominalValues();
83 MEB::OutArgs<Scalar> outArgs = this->appModel_->createOutArgs();
84 inArgs.set_x(stageX);
85 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
86 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
87 outArgs.set_f(stageXDot[0]); // K1
88 this->appModel_->evalModel(inArgs, outArgs);
89
90 this->stepperErrorNormCalculator_->setRelativeTolerance(errorRel);
91 this->stepperErrorNormCalculator_->setAbsoluteTolerance(errorAbs);
92
93 Scalar d0 = this->stepperErrorNormCalculator_->errorNorm(stageX);
94 Scalar d1 = this->stepperErrorNormCalculator_->errorNorm(stageXDot[0]);
95
96 // b) first guess for the step size
97 dt = Teuchos::as<Scalar>(0.01)*(d0/d1);
98
99 // c) perform one explicit Euler step (X_1)
100 Thyra::Vp_StV(stageX.ptr(), dt, *(stageXDot[0]));
101
102 // compute F(t_0 + dt, X_1)
103 inArgs.set_x(stageX);
104 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time + dt);
105 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
106 outArgs.set_f(stageXDot[1]); // K2
107 this->appModel_->evalModel(inArgs, outArgs);
108
109 // d) compute estimate of the second derivative of the solution
110 // d2 = || f(t_0 + dt, X_1) - f(t_0, X_0) || / dt
111 Teuchos::RCP<Thyra::VectorBase<Scalar> > errX;
112 errX = Thyra::createMember(this->appModel_->get_f_space());
113 assign(errX.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
114 Thyra::V_VmV(errX.ptr(), *(stageXDot[1]), *(stageXDot[0]));
115 Scalar d2 = this->stepperErrorNormCalculator_->errorNorm(errX) / dt;
116
117 // e) compute step size h_1 (from m = 0 order Taylor series)
118 Scalar max_d1_d2 = std::max(d1, d2);
119 Scalar h1 = std::pow((0.01/max_d1_d2),(1.0/(order+1)));
120
121 // f) propose starting step size
122 dt = std::min(100*dt, h1);
123 return dt;
124}
125
126
127template<class Scalar>
128Teuchos::RCP<const Teuchos::ParameterList>
130{
131 return this->getValidParametersBasicERK();
132}
133
134
135template<class Scalar>
136Teuchos::RCP<Teuchos::ParameterList>
138{
139 auto pl = this->getValidParametersBasic();
140 pl->template set<bool>("Use Embedded", this->getUseEmbedded(),
141 "'Whether to use Embedded Stepper (if available) or not\n"
142 " 'true' - Stepper will compute embedded solution and is adaptive.\n"
143 " 'false' - Stepper is not embedded(adaptive).\n");
144 pl->template set<std::string>("Description", this->getDescription());
145
146 return pl;
147}
148
149
150template<class Scalar>
152{
153 TEUCHOS_TEST_FOR_EXCEPTION( this->tableau_ == Teuchos::null, std::logic_error,
154 "Error - Need to set the tableau, before calling "
155 "StepperExplicitRK::initialize()\n");
156
157 TEUCHOS_TEST_FOR_EXCEPTION( this->appModel_==Teuchos::null, std::logic_error,
158 "Error - Need to set the model, setModel(), before calling "
159 "StepperExplicitRK::initialize()\n");
160
162}
163
164
165template<class Scalar>
167 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
168{
170
171 // Set the stage vectors
172 int numStages = this->tableau_->numStages();
173 stageXDot_.resize(numStages);
174 for (int i=0; i<numStages; ++i) {
175 stageXDot_[i] = Thyra::createMember(this->appModel_->get_f_space());
176 assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
177 }
178
179 this->setEmbeddedMemory();
180 this->setErrorNorm();
181
182 this->isInitialized_ = false;
183}
184
185
186template<class Scalar>
188{
189 if (this->getModel() == Teuchos::null)
190 return; // Embedded memory will be set when setModel() is called.
191
192 if ( this->tableau_->isEmbedded() && this->getUseEmbedded() ){
193 this->ee_ = Thyra::createMember(this->appModel_->get_f_space());
194 this->abs_u0 = Thyra::createMember(this->appModel_->get_f_space());
195 this->abs_u = Thyra::createMember(this->appModel_->get_f_space());
196 this->sc = Thyra::createMember(this->appModel_->get_f_space());
197 } else {
198 this->ee_ = Teuchos::null;
199 this->abs_u0 = Teuchos::null;
200 this->abs_u = Teuchos::null;
201 this->sc = Teuchos::null;
202 }
203}
204
205
206template<class Scalar>
208 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
209{
210 if (this->getUseFSAL())
211 this->setStepperXDot(stageXDot_.back());
212 else
213 this->setStepperXDot(stageXDot_[0]);
214
216
217 auto xDot = solutionHistory->getCurrentState()->getXDot();
218 if (xDot != Teuchos::null && this->getUseFSAL())
219 Thyra::assign(this->getStepperXDot().ptr(), *(xDot));
220}
221
222
223template<class Scalar>
225 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
226{
227 this->checkInitialized();
228
229 using Teuchos::RCP;
230
231 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperExplicitRK::takeStep()");
232 {
233 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
234 std::logic_error,
235 "Error - StepperExplicitRK<Scalar>::takeStep(...)\n"
236 "Need at least two SolutionStates for ExplicitRK.\n"
237 " Number of States = " << solutionHistory->getNumStates() << "\n"
238 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
239 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
240
241 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
242 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
243 const Scalar dt = workingState->getTimeStep();
244 const Scalar time = currentState->getTime();
245
246 const int numStages = this->tableau_->numStages();
247 Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
248 Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
249 Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
250
251 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
252
253 RCP<StepperExplicitRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
254 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
256
257 // Compute stage solutions
258 for (int i=0; i < numStages; ++i) {
259 this->setStageNumber(i);
260 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
261 for (int j=0; j < i; ++j) {
262 if (A(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
263 Thyra::Vp_StV(workingState->getX().ptr(), dt*A(i,j), *stageXDot_[j]);
264 }
265 }
266 this->setStepperXDot(stageXDot_[i]);
267
268 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
270 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
272 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
274 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
276
277 if ( i == 0 && this->getUseFSAL() &&
278 workingState->getNConsecutiveFailures() == 0 ) {
279 RCP<Thyra::VectorBase<Scalar> > tmp = stageXDot_[0];
280 stageXDot_[0] = stageXDot_.back();
281 stageXDot_.back() = tmp;
282 this->setStepperXDot(stageXDot_[0]);
283
284 } else {
285 const Scalar ts = time + c(i)*dt;
286 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(dt));
287
288 // Evaluate xDot = f(x,t).
289 this->evaluateExplicitODE(stageXDot_[i], workingState->getX(), ts, p);
290 }
291
292 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
294 }
295 // reset the stage number
296 this->setStageNumber(-1);
297
298 // Sum for solution: x_n = x_n-1 + Sum{ b(i) * dt*f(i) }
299 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
300 for (int i=0; i < numStages; ++i) {
301 if (b(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
302 Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
303 }
304 }
305
306 if (this->getUseFSAL()) {
307 if (numStages == 1) {
308 const Scalar ts = time + dt;
309 auto p = Teuchos::rcp(new ExplicitODEParameters<Scalar>(dt));
310 // Evaluate xDot = f(x,t).
311 this->evaluateExplicitODE(stageXDot_[0], workingState->getX(), ts, p);
312 }
313 if (workingState->getXDot() != Teuchos::null)
314 Thyra::assign((workingState->getXDot()).ptr(), *(stageXDot_.back()));
315 }
316
317
318 // At this point, the stepper has passed.
319 // But when using adaptive time stepping, the embedded method
320 // can change the step status
321 workingState->setSolutionStatus(Status::PASSED);
322
323 if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
324
325 const Scalar tolRel = workingState->getTolRel();
326 const Scalar tolAbs = workingState->getTolAbs();
327
328 // update the tolerance
329 this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
330 this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
331
332 // just compute the error weight vector
333 // (all that is needed is the error, and not the embedded solution)
334 Teuchos::SerialDenseVector<int,Scalar> errWght = b ;
335 errWght -= this->tableau_->bstar();
336
337 // Compute local truncation error estimate: | u^{n+1} - \hat{u}^{n+1} |
338 // Sum for solution: ee_n = Sum{ (b(i) - bstar(i)) * dt*f(i) }
339 assign(this->ee_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
340 for (int i=0; i < numStages; ++i) {
341 if (errWght(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
342 Thyra::Vp_StV(this->ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
343 }
344 }
345
346
347 Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
348 workingState->setErrorRel(err);
349
350 // Test if step should be rejected
351 if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
352 workingState->setSolutionStatus(Status::FAILED);
353 }
354
355 workingState->setOrder(this->getOrder());
356 workingState->computeNorms(currentState);
357 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
359 }
360 return;
361}
362
363
370template<class Scalar>
371Teuchos::RCP<Tempus::StepperState<Scalar> > StepperExplicitRK<Scalar>::
373{
374 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
375 rcp(new StepperState<Scalar>(this->getStepperType()));
376 return stepperState;
377}
378
379
380template<class Scalar>
382 Teuchos::FancyOStream &out,
383 const Teuchos::EVerbosityLevel verbLevel) const
384{
385 out.setOutputToRootOnly(0);
386 out << std::endl;
387 Stepper<Scalar>::describe(out, verbLevel);
388 StepperExplicit<Scalar>::describe(out, verbLevel);
389
390 out << "--- StepperExplicitRK ---\n";
391 if (this->tableau_ != Teuchos::null) this->tableau_->describe(out, verbLevel);
392 out << " tableau_ = " << this->tableau_ << std::endl;
393 out << " stepperRKAppAction_= " << this->stepperRKAppAction_ << std::endl;
394 out << " stageXDot_.size() = " << stageXDot_.size() << std::endl;
395 const int numStages = stageXDot_.size();
396 for (int i=0; i<numStages; ++i)
397 out << " stageXDot_["<<i<<"] = " << stageXDot_[i] << std::endl;
398 out << " useEmbedded_ = "
399 << Teuchos::toString(this->useEmbedded_) << std::endl;
400 out << " ee_ = " << this->ee_ << std::endl;
401 out << " abs_u0 = " << this->abs_u0 << std::endl;
402 out << " abs_u = " << this->abs_u << std::endl;
403 out << " sc = " << this->sc << std::endl;
404 out << "-------------------------" << std::endl;
405}
406
407
408template<class Scalar>
409bool StepperExplicitRK<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
410{
411 out.setOutputToRootOnly(0);
412 bool isValidSetup = true;
413
414 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
415 if ( !StepperExplicit<Scalar>::isValidSetup(out) ) isValidSetup = false;
416
417 if (this->tableau_ == Teuchos::null) {
418 isValidSetup = false;
419 out << "The tableau is not set!\n";
420 }
421
422 if (this->stepperRKAppAction_ == Teuchos::null) {
423 isValidSetup = false;
424 out << "The AppAction is not set!\n";
425 }
426
427 return isValidSetup;
428}
429
430
431} // namespace Tempus
432#endif // Tempus_StepperExplicitRK_impl_hpp
virtual void setup(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel, bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck, bool useEmbedded, const Teuchos::RCP< StepperRKAppAction< Scalar > > &stepperRKAppAction)
Setup for constructor.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
virtual Scalar getInitTimeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory) const
virtual void setupDefault()
Default setup for constructor.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set model.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Get a default (initial) StepperState.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< Teuchos::ParameterList > getValidParametersBasicERK() const
virtual void initialize()
Initialize during construction and after changing input parameters.
Thyra Base interface for implicit time steppers.
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.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Application Action for StepperRKBase.
Thyra Base interface for time steppers.
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const