148 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
150 this->checkInitialized();
154 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperDIRK::takeStep()");
156 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
158 "Error - StepperDIRK<Scalar>::takeStep(...)\n"
159 "Need at least two SolutionStates for DIRK.\n"
160 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
161 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
162 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
164 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
165 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
166 const Scalar dt = workingState->getTimeStep();
167 const Scalar time = currentState->getTime();
169 const int numStages = this->tableau_->numStages();
170 Teuchos::SerialDenseMatrix<int,Scalar> A = this->tableau_->A();
171 Teuchos::SerialDenseVector<int,Scalar> b = this->tableau_->b();
172 Teuchos::SerialDenseVector<int,Scalar> c = this->tableau_->c();
175 if ( this->getResetInitialGuess() && (!this->getZeroInitialGuess()) )
176 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
178 RCP<StepperDIRK<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
179 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
184 Thyra::SolveStatus<Scalar> sStatus;
185 for (
int i=0; i < numStages; ++i) {
186 this->setStageNumber(i);
188 Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
189 for (
int j=0; j < i; ++j) {
190 if (A(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
191 Thyra::Vp_StV(xTilde_.ptr(), dt*A(i,j), *(stageXDot_[j]));
194 this->setStepperXDot(stageXDot_[i]);
196 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
199 Scalar ts = time + c(i)*dt;
202 bool isNeeded =
false;
203 for (
int k=i+1; k<numStages; ++k)
if (A(k,i) != 0.0) isNeeded =
true;
204 if (b(i) != 0.0) isNeeded =
true;
205 if (this->tableau_->isEmbedded() && this->getUseEmbedded() &&
206 this->tableau_->bstar()(i) != 0.0)
208 if (isNeeded ==
false) {
209 assign(stageXDot_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
210 }
else if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
212 if (i == 0 && this->getUseFSAL() &&
213 workingState->getNConsecutiveFailures() == 0) {
215 RCP<Thyra::VectorBase<Scalar> > tmp = stageXDot_[0];
216 stageXDot_[0] = stageXDot_.back();
217 stageXDot_.back() = tmp;
218 this->setStepperXDot(stageXDot_[0]);
221 typedef Thyra::ModelEvaluatorBase MEB;
222 MEB::InArgs<Scalar> inArgs = this->wrapperModel_->getInArgs();
223 MEB::OutArgs<Scalar> outArgs = this->wrapperModel_->getOutArgs();
224 inArgs.set_x(xTilde_);
225 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(ts);
226 if (inArgs.supports(MEB::IN_ARG_x_dot))
227 inArgs.set_x_dot(Teuchos::null);
228 outArgs.set_f(stageXDot_[i]);
230 this->wrapperModel_->getAppModel()->evalModel(inArgs,outArgs);
234 const Scalar alpha = 1.0/(dt*A(i,i));
235 const Scalar beta = 1.0;
238 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
239 Teuchos::rcp(
new StepperDIRKTimeDerivative<Scalar>(
240 alpha,xTilde_.getConst()));
242 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
245 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
248 sStatus = this->solveImplicitODE(workingState->getX(), stageXDot_[i], ts, p);
250 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass=
false;
252 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
255 timeDer->compute(workingState->getX(), stageXDot_[i]);
257 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
259 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
263 this->setStageNumber(-1);
266 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
267 for (
int i=0; i < numStages; ++i) {
268 if (b(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
269 Thyra::Vp_StV((workingState->getX()).ptr(), dt*b(i), *(stageXDot_[i]));
273 if (this->tableau_->isEmbedded() && this->getUseEmbedded()) {
274 const Scalar tolRel = workingState->getTolRel();
275 const Scalar tolAbs = workingState->getTolAbs();
278 this->stepperErrorNormCalculator_->setRelativeTolerance(tolRel);
279 this->stepperErrorNormCalculator_->setAbsoluteTolerance(tolAbs);
283 Teuchos::SerialDenseVector<int,Scalar> errWght = b ;
284 errWght -= this->tableau_->bstar();
288 assign(this->ee_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
289 for (
int i=0; i < numStages; ++i) {
290 if (errWght(i) != Teuchos::ScalarTraits<Scalar>::zero()) {
291 Thyra::Vp_StV(this->ee_.ptr(), dt*errWght(i), *(stageXDot_[i]));
295 Scalar err = this->stepperErrorNormCalculator_->computeWRMSNorm(currentState->getX(), workingState->getX(), this->ee_);
296 workingState->setErrorRel(err);
299 if (std::isinf(err) || std::isnan(err) || err > Teuchos::as<Scalar>(1.0))
306 workingState->setOrder(this->getOrder());
307 workingState->computeNorms(currentState);
308 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,