107 std::string stepperType,
111 if (stepperType ==
"") stepperType =
"Partitioned IMEX RK SSP2";
113 if (stepperType ==
"Partitioned IMEX RK 1st order") {
116 typedef Teuchos::ScalarTraits<Scalar> ST;
118 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
119 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
120 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
121 const Scalar one = ST::one();
122 const Scalar zero = ST::zero();
141 "Explicit Tableau - Partitioned IMEX RK 1st order", A, b, c, order,
143 expTableau->setTVD(
true);
144 expTableau->setTVDCoeff(2.0);
146 this->setExplicitTableau(expTableau);
150 typedef Teuchos::ScalarTraits<Scalar> ST;
152 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
153 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
154 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
155 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
156 const Scalar one = ST::one();
157 const Scalar zero = ST::zero();
176 "Implicit Tableau - Partitioned IMEX RK 1st order", A, b, c, order,
178 impTableau->setTVD(
true);
179 impTableau->setTVDCoeff(sspcoef);
181 this->setImplicitTableau(impTableau);
183 this->setStepperName(
"Partitioned IMEX RK 1st order");
184 this->setStepperType(
"Partitioned IMEX RK 1st order");
187 else if (stepperType ==
"Partitioned IMEX RK SSP2") {
190 this->setExplicitTableau(stepperERK->getTableau());
194 stepperSDIRK->setGammaType(
"2nd Order L-stable");
195 this->setImplicitTableau(stepperSDIRK->getTableau());
197 this->setStepperName(
"Partitioned IMEX RK SSP2");
198 this->setStepperType(
"Partitioned IMEX RK SSP2");
201 else if (stepperType ==
"Partitioned IMEX RK ARS 233") {
202 typedef Teuchos::ScalarTraits<Scalar> ST;
204 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
205 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
206 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
207 const Scalar one = ST::one();
208 const Scalar zero = ST::zero();
209 const Scalar onehalf = ST::one() / (2 * ST::one());
210 const Scalar gamma = (3 * one + ST::squareroot(3 * one)) / (6 * one);
220 A(2, 0) = (gamma - 1.0);
221 A(2, 1) = (2.0 - 2.0 * gamma);
237 "Explicit Tableau - Partitioned IMEX RK ARS 233", A, b, c, order,
240 this->setExplicitTableau(expTableau);
252 A(2, 1) = (1.0 - 2.0 * gamma);
268 "Implicit Tableau - Partitioned IMEX RK ARS 233", A, b, c, order,
271 this->setImplicitTableau(impTableau);
273 this->setStepperName(
"Partitioned IMEX RK ARS 233");
274 this->setStepperType(
"Partitioned IMEX RK ARS 233");
277 else if (stepperType ==
"General Partitioned IMEX RK") {
278 if (explicitTableau == Teuchos::null) {
281 this->setExplicitTableau(stepperERK->getTableau());
284 this->setExplicitTableau(explicitTableau);
287 if (implicitTableau == Teuchos::null) {
291 stepperSDIRK->setGammaType(
"2nd Order L-stable");
292 this->setImplicitTableau(stepperSDIRK->getTableau());
295 this->setImplicitTableau(implicitTableau);
298 this->setStepperName(
"General Partitioned IMEX RK");
299 this->setStepperType(
"General Partitioned IMEX RK");
303 TEUCHOS_TEST_FOR_EXCEPTION(
304 true, std::logic_error,
305 "Error - Not a valid StepperIMEX_RK_Partition type! Stepper Type = "
306 << stepperType <<
"\n"
307 <<
" Current valid types are: \n"
308 <<
" 'Partitioned IMEX RK 1st order'\n"
309 <<
" 'Partitioned IMEX RK SSP2'\n"
310 <<
" 'Partitioned IMEX RK ARS 233'\n"
311 <<
" 'General Partitioned IMEX RK'\n");
314 TEUCHOS_TEST_FOR_EXCEPTION(
315 explicitTableau_ == Teuchos::null, std::runtime_error,
316 "Error - StepperIMEX_RK_Partition - Explicit tableau is null!");
317 TEUCHOS_TEST_FOR_EXCEPTION(
318 implicitTableau_ == Teuchos::null, std::runtime_error,
319 "Error - StepperIMEX_RK_Partition - Implicit tableau is null!");
320 TEUCHOS_TEST_FOR_EXCEPTION(
321 explicitTableau_->numStages() != implicitTableau_->numStages(),
323 "Error - StepperIMEX_RK_Partition - Number of stages do not match!\n"
324 <<
" Explicit tableau = " << explicitTableau_->description() <<
"\n"
325 <<
" number of stages = " << explicitTableau_->numStages() <<
"\n"
326 <<
" Implicit tableau = " << implicitTableau_->description() <<
"\n"
327 <<
" number of stages = " << implicitTableau_->numStages()
330 this->isInitialized_ =
false;
335 Teuchos::RCP<Teuchos::ParameterList> pl, std::string stepperType)
338 if (stepperType ==
"") {
339 if (pl == Teuchos::null)
340 stepperType =
"Partitioned IMEX RK SSP2";
343 pl->get<std::string>(
"Stepper Type",
"Partitioned IMEX RK SSP2");
346 if (stepperType !=
"General Partitioned IMEX RK") {
347 this->setTableaus(stepperType);
350 if (pl != Teuchos::null) {
351 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau;
352 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau;
353 if (pl->isSublist(
"IMEX-RK Explicit Stepper")) {
354 RCP<Teuchos::ParameterList> explicitPL =
355 Teuchos::rcp(
new Teuchos::ParameterList(
356 pl->sublist(
"IMEX-RK Explicit Stepper")));
358 auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
359 auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> >(
361 TEUCHOS_TEST_FOR_EXCEPTION(
362 stepperERK == Teuchos::null, std::logic_error,
363 "Error - The explicit component of a general partitioned IMEX RK "
364 "stepper was not specified as an ExplicitRK stepper");
365 explicitTableau = stepperERK->getTableau();
368 if (pl->isSublist(
"IMEX-RK Implicit Stepper")) {
369 RCP<Teuchos::ParameterList> implicitPL =
370 Teuchos::rcp(
new Teuchos::ParameterList(
371 pl->sublist(
"IMEX-RK Implicit Stepper")));
373 auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
375 Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> >(stepperTemp,
true);
376 TEUCHOS_TEST_FOR_EXCEPTION(
377 stepperDIRK == Teuchos::null, std::logic_error,
378 "Error - The implicit component of a general partitioned IMEX RK "
379 "stepper was not specified as an DIRK stepper");
380 implicitTableau = stepperDIRK->getTableau();
383 TEUCHOS_TEST_FOR_EXCEPTION(
384 !(explicitTableau != Teuchos::null &&
385 implicitTableau != Teuchos::null),
387 "Error - A parameter list was used to setup a general partitioned "
388 "IMEX RK stepper, but did not "
389 "specify both an explicit and an implicit tableau!\n");
391 this->setTableaus(stepperType, explicitTableau, implicitTableau);
393 this->setOrder(pl->get<
int>(
"overall order", 1));
648 this->checkInitialized();
651 using Teuchos::SerialDenseMatrix;
652 using Teuchos::SerialDenseVector;
654 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperIMEX_RK_Partition::takeStep()");
656 TEUCHOS_TEST_FOR_EXCEPTION(
657 solutionHistory->getNumStates() < 2, std::logic_error,
658 "Error - StepperIMEX_RK_Partition<Scalar>::takeStep(...)\n"
659 <<
"Need at least two SolutionStates for IMEX_RK_Partition.\n"
660 <<
" Number of States = " << solutionHistory->getNumStates()
661 <<
"\nTry setting in \"Solution History\" \"Storage Type\" = "
662 <<
"\"Undo\"\n or \"Storage Type\" = \"Static\" and \"Storage Limit\" = "
665 RCP<SolutionState<Scalar> > currentState =
666 solutionHistory->getCurrentState();
667 RCP<SolutionState<Scalar> > workingState =
668 solutionHistory->getWorkingState();
669 const Scalar dt = workingState->getTimeStep();
670 const Scalar time = currentState->getTime();
672 const int numStages = explicitTableau_->numStages();
673 const SerialDenseMatrix<int, Scalar>& AHat = explicitTableau_->A();
674 const SerialDenseVector<int, Scalar>& bHat = explicitTableau_->b();
675 const SerialDenseVector<int, Scalar>& cHat = explicitTableau_->c();
676 const SerialDenseMatrix<int, Scalar>& A = implicitTableau_->A();
677 const SerialDenseVector<int, Scalar>& b = implicitTableau_->b();
678 const SerialDenseVector<int, Scalar>& c = implicitTableau_->c();
680 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
681 wrapperModelPairPartIMEX = Teuchos::rcp_dynamic_cast<
683 this->wrapperModel_);
686 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
687 RCP<Thyra::VectorBase<Scalar> > stageY =
688 wrapperModelPairPartIMEX->getExplicitOnlyVector(workingState->getX());
689 RCP<Thyra::VectorBase<Scalar> > stageX =
690 wrapperModelPairPartIMEX->getIMEXVector(workingState->getX());
692 RCP<StepperIMEX_RK_Partition<Scalar> > thisStepper =
693 Teuchos::rcpFromRef(*
this);
694 this->stepperRKAppAction_->execute(
695 solutionHistory, thisStepper,
699 for (
int i = 0; i < numStages; ++i) {
700 this->setStageNumber(i);
702 Thyra::assign(stageY.ptr(),
703 *(wrapperModelPairPartIMEX->getExplicitOnlyVector(
704 currentState->getX())));
707 *(wrapperModelPairPartIMEX->getIMEXVector(currentState->getX())));
708 for (
int j = 0; j < i; ++j) {
709 if (AHat(i, j) != Teuchos::ScalarTraits<Scalar>::zero()) {
710 RCP<Thyra::VectorBase<Scalar> > stageFy =
711 wrapperModelPairPartIMEX->getExplicitOnlyVector(stageF_[j]);
712 RCP<Thyra::VectorBase<Scalar> > stageFx =
713 wrapperModelPairPartIMEX->getIMEXVector(stageF_[j]);
714 Thyra::Vp_StV(stageY.ptr(), -dt * AHat(i, j), *stageFy);
715 Thyra::Vp_StV(xTilde_.ptr(), -dt * AHat(i, j), *stageFx);
717 if (A(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
718 Thyra::Vp_StV(xTilde_.ptr(), -dt * A(i, j), *(stageGx_[j]));
721 this->stepperRKAppAction_->execute(
722 solutionHistory, thisStepper,
725 Scalar ts = time + c(i) * dt;
726 Scalar tHats = time + cHat(i) * dt;
727 if (A(i, i) == Teuchos::ScalarTraits<Scalar>::zero()) {
729 bool isNeeded =
false;
730 for (
int k = i + 1; k < numStages; ++k)
731 if (A(k, i) != 0.0) isNeeded =
true;
732 if (b(i) != 0.0) isNeeded =
true;
733 if (isNeeded ==
false) {
735 assign(stageGx_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
738 Thyra::assign(stageX.ptr(), *xTilde_);
739 evalImplicitModelExplicitly(stageX, stageY, ts, dt, i, stageGx_[i]);
744 const Scalar alpha = Scalar(1.0) / (dt * A(i, i));
745 const Scalar beta = Scalar(1.0);
748 RCP<TimeDerivative<Scalar> > timeDer =
750 alpha, xTilde_.getConst()));
755 this->stepperRKAppAction_->execute(
756 solutionHistory, thisStepper,
759 wrapperModelPairPartIMEX->setUseImplicitModel(
true);
760 this->solver_->setModel(wrapperModelPairPartIMEX);
762 Thyra::SolveStatus<Scalar> sStatus;
763 if (wrapperModelPairPartIMEX->getParameterIndex() >= 0)
764 sStatus = this->solveImplicitODE(
765 stageX, stageGx_[i], ts, p, stageY,
766 wrapperModelPairPartIMEX->getParameterIndex());
768 sStatus = this->solveImplicitODE(stageX, stageGx_[i], ts, p);
770 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass =
false;
772 wrapperModelPairPartIMEX->setUseImplicitModel(
false);
774 this->stepperRKAppAction_->execute(
775 solutionHistory, thisStepper,
779 Thyra::V_StVpStV(stageGx_[i].ptr(), -alpha, *stageX, alpha, *xTilde_);
782 this->stepperRKAppAction_->execute(
783 solutionHistory, thisStepper,
785 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
786 this->stepperRKAppAction_->execute(
787 solutionHistory, thisStepper,
793 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
794 RCP<Thyra::VectorBase<Scalar> > Z = workingState->getX();
795 RCP<Thyra::VectorBase<Scalar> > X =
796 wrapperModelPairPartIMEX->getIMEXVector(Z);
797 for (
int i = 0; i < numStages; ++i) {
798 if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
799 Thyra::Vp_StV(Z.ptr(), -dt * bHat(i), *(stageF_[i]));
800 if (b(i) != Teuchos::ScalarTraits<Scalar>::zero())
801 Thyra::Vp_StV(X.ptr(), -dt * b(i), *(stageGx_[i]));
808 workingState->setOrder(this->getOrder());
809 workingState->computeNorms(currentState);
810 this->stepperRKAppAction_->execute(
811 solutionHistory, thisStepper,
815 this->setStageNumber(-1);