57 const Teuchos::RCP<SolutionHistory<Scalar> >& sh)
const
60 Scalar dt = Scalar(1.0e+99);
61 if (!this->getUseEmbedded())
return dt;
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();
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()));
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());
81 typedef Thyra::ModelEvaluatorBase MEB;
82 MEB::InArgs<Scalar> inArgs = this->appModel_->getNominalValues();
83 MEB::OutArgs<Scalar> outArgs = this->appModel_->createOutArgs();
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]);
88 this->appModel_->evalModel(inArgs, outArgs);
90 this->stepperErrorNormCalculator_->setRelativeTolerance(errorRel);
91 this->stepperErrorNormCalculator_->setAbsoluteTolerance(errorAbs);
93 Scalar d0 = this->stepperErrorNormCalculator_->errorNorm(stageX);
94 Scalar d1 = this->stepperErrorNormCalculator_->errorNorm(stageXDot[0]);
97 dt = Teuchos::as<Scalar>(0.01)*(d0/d1);
100 Thyra::Vp_StV(stageX.ptr(), dt, *(stageXDot[0]));
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]);
107 this->appModel_->evalModel(inArgs, outArgs);
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;
118 Scalar max_d1_d2 = std::max(d1, d2);
119 Scalar h1 = std::pow((0.01/max_d1_d2),(1.0/(order+1)));
122 dt = std::min(100*dt, h1);
225 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
227 this->checkInitialized();
231 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperExplicitRK::takeStep()");
233 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
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");
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();
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();
251 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
253 RCP<StepperExplicitRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
254 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
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]);
266 this->setStepperXDot(stageXDot_[i]);
268 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
270 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
272 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
274 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
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]);
285 const Scalar ts = time + c(i)*dt;
286 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(dt));
289 this->evaluateExplicitODE(stageXDot_[i], workingState->getX(), ts, p);
292 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
296 this->setStageNumber(-1);
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]));
306 if (this->getUseFSAL()) {
307 if (numStages == 1) {
308 const Scalar ts = time + dt;
309 auto p = Teuchos::rcp(
new ExplicitODEParameters<Scalar>(dt));
311 this->evaluateExplicitODE(stageXDot_[0], workingState->getX(), ts, p);
313 if (workingState->getXDot() != Teuchos::null)
314 Thyra::assign((workingState->getXDot()).ptr(), *(stageXDot_.back()));
323 if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
325 const Scalar tolRel = workingState->getTolRel();
326 const Scalar tolAbs = workingState->getTolAbs();
329 this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
330 this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
334 Teuchos::SerialDenseVector<int,Scalar> errWght = b ;
335 errWght -= this->tableau_->bstar();
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]));
347 Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
348 workingState->setErrorRel(err);
351 if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
355 workingState->setOrder(this->getOrder());
356 workingState->computeNorms(currentState);
357 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
382 Teuchos::FancyOStream &out,
383 const Teuchos::EVerbosityLevel verbLevel)
const
385 out.setOutputToRootOnly(0);
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;