124 std::string stepperType,
128 if (stepperType ==
"") stepperType =
"IMEX RK SSP2";
130 if (stepperType ==
"IMEX RK 1st order") {
133 typedef Teuchos::ScalarTraits<Scalar> ST;
135 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
136 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
137 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
138 const Scalar one = ST::one();
139 const Scalar zero = ST::zero();
157 auto expTableau = Teuchos::rcp(
159 A, b, c, order, order, order));
160 expTableau->setTVD(
true);
161 expTableau->setTVDCoeff(2.0);
163 this->setExplicitTableau(expTableau);
167 typedef Teuchos::ScalarTraits<Scalar> ST;
169 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
170 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
171 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
172 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
173 const Scalar one = ST::one();
174 const Scalar zero = ST::zero();
192 auto impTableau = Teuchos::rcp(
194 A, b, c, order, order, order));
195 impTableau->setTVD(
true);
196 impTableau->setTVDCoeff(sspcoef);
198 this->setImplicitTableau(impTableau);
200 this->setStepperName(
"IMEX RK 1st order");
201 this->setStepperType(
"IMEX RK 1st order");
204 else if (stepperType ==
"SSP1_111") {
207 typedef Teuchos::ScalarTraits<Scalar> ST;
208 const int NumStages = 1;
210 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
211 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
212 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
213 const Scalar one = ST::one();
214 const Scalar zero = ST::zero();
226 "Explicit Tableau - SSP1_111", A, b, c, order, order, order));
227 expTableau->setTVD(
true);
228 expTableau->setTVDCoeff(1.0);
230 this->setExplicitTableau(expTableau);
234 typedef Teuchos::ScalarTraits<Scalar> ST;
235 const int NumStages = 1;
237 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
238 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
239 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
240 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
241 const Scalar one = ST::one();
253 "Implicit Tableau - SSP1_111", A, b, c, order, order, order));
254 impTableau->setTVD(
true);
255 impTableau->setTVDCoeff(sspcoef);
257 this->setImplicitTableau(impTableau);
259 this->setStepperName(
"SSP1_111");
260 this->setStepperType(
"SSP1_111");
263 else if (stepperType ==
"IMEX RK SSP2" || stepperType ==
"SSP2_222_L") {
266 this->setExplicitTableau(stepperERK->getTableau());
270 stepperSDIRK->setGammaType(
"2nd Order L-stable");
271 this->setImplicitTableau(stepperSDIRK->getTableau());
273 this->setStepperName(
"IMEX RK SSP2");
274 this->setStepperType(
"IMEX RK SSP2");
277 else if (stepperType ==
"SSP2_222" || stepperType ==
"SSP2_222_A") {
280 this->setExplicitTableau(stepperERK->getTableau());
284 stepperSDIRK->setGammaType(
"gamma");
285 stepperSDIRK->setGamma(0.5);
286 this->setImplicitTableau(stepperSDIRK->getTableau());
288 this->setStepperName(
"SSP2_222");
289 this->setStepperType(
"SSP2_222");
292 else if (stepperType ==
"IMEX RK SSP3" || stepperType ==
"SSP3_332") {
295 this->setExplicitTableau(stepperERK->getTableau());
299 this->setImplicitTableau(stepperSDIRK->getTableau());
301 this->setStepperName(
"IMEX RK SSP3");
302 this->setStepperType(
"IMEX RK SSP3");
305 else if (stepperType ==
"IMEX RK ARS 233" || stepperType ==
"ARS 233") {
306 typedef Teuchos::ScalarTraits<Scalar> ST;
308 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
309 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
310 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
311 const Scalar one = ST::one();
312 const Scalar zero = ST::zero();
313 const Scalar onehalf = ST::one() / (2 * ST::one());
314 const Scalar gamma = (3 * one + ST::squareroot(3 * one)) / (6 * one);
324 A(2, 0) = (gamma - 1.0);
325 A(2, 1) = (2.0 - 2.0 * gamma);
341 "Partition IMEX-RK Explicit Stepper", A, b, c, order, order, order));
343 this->setExplicitTableau(expTableau);
355 A(2, 1) = (1.0 - 2.0 * gamma);
371 "Partition IMEX-RK Implicit Stepper", A, b, c, order, order, order));
373 this->setImplicitTableau(impTableau);
375 this->setStepperName(
"IMEX RK ARS 233");
376 this->setStepperType(
"IMEX RK ARS 233");
379 else if (stepperType ==
"General IMEX RK") {
380 if (explicitTableau == Teuchos::null) {
383 this->setExplicitTableau(stepperERK->getTableau());
386 this->setExplicitTableau(explicitTableau);
389 if (explicitTableau == Teuchos::null) {
393 stepperSDIRK->setGammaType(
"2nd Order L-stable");
394 this->setImplicitTableau(stepperSDIRK->getTableau());
397 this->setImplicitTableau(implicitTableau);
400 this->setStepperName(
"General IMEX RK");
401 this->setStepperType(
"General IMEX RK");
405 TEUCHOS_TEST_FOR_EXCEPTION(
406 true, std::logic_error,
407 "Error - Not a valid StepperIMEX_RK type! Stepper Type = "
408 << stepperType <<
"\n"
409 <<
" Current valid types are: \n"
410 <<
" 'IMEX RK 1st order\n"
412 <<
" 'IMEX RK SSP2' ('SSP2_222_L')\n"
413 <<
" 'SSP2_222' ('SSP2_222_A')\n"
414 <<
" 'IMEX RK SSP3' ('SSP3_332')\n"
415 <<
" 'IMEX RK ARS 233' ('ARS 233')\n"
416 <<
" 'General IMEX RK'\n");
419 TEUCHOS_TEST_FOR_EXCEPTION(
420 explicitTableau_ == Teuchos::null, std::runtime_error,
421 "Error - StepperIMEX_RK - Explicit tableau is null!");
422 TEUCHOS_TEST_FOR_EXCEPTION(
423 implicitTableau_ == Teuchos::null, std::runtime_error,
424 "Error - StepperIMEX_RK - Implicit tableau is null!");
425 TEUCHOS_TEST_FOR_EXCEPTION(
426 explicitTableau_->numStages() != implicitTableau_->numStages(),
428 "Error - StepperIMEX_RK - Number of stages do not match!\n"
429 <<
" Explicit tableau = " << explicitTableau_->description() <<
"\n"
430 <<
" number of stages = " << explicitTableau_->numStages() <<
"\n"
431 <<
" Implicit tableau = " << implicitTableau_->description() <<
"\n"
432 <<
" number of stages = " << implicitTableau_->numStages() <<
"\n");
434 this->isInitialized_ =
false;
439 Teuchos::RCP<Teuchos::ParameterList> pl, std::string stepperType)
442 if (stepperType ==
"") {
443 if (pl == Teuchos::null)
444 stepperType =
"IMEX RK SSP2";
446 stepperType = pl->get<std::string>(
"Stepper Type",
"IMEX RK SSP2");
449 if (stepperType !=
"General IMEX RK") {
450 this->setTableaus(stepperType);
453 if (pl != Teuchos::null) {
454 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau;
455 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau;
456 if (pl->isSublist(
"IMEX-RK Explicit Stepper")) {
457 RCP<Teuchos::ParameterList> explicitPL =
458 Teuchos::rcp(
new Teuchos::ParameterList(
459 pl->sublist(
"IMEX-RK Explicit Stepper")));
461 auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
462 auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> >(
464 TEUCHOS_TEST_FOR_EXCEPTION(
465 stepperERK == Teuchos::null, std::logic_error,
466 "Error - The explicit component of a general IMEX RK stepper was "
467 "not specified as an ExplicitRK stepper");
468 explicitTableau = stepperERK->getTableau();
471 if (pl->isSublist(
"IMEX-RK Implicit Stepper")) {
472 RCP<Teuchos::ParameterList> implicitPL =
473 Teuchos::rcp(
new Teuchos::ParameterList(
474 pl->sublist(
"IMEX-RK Implicit Stepper")));
476 auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
478 Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> >(stepperTemp,
true);
479 TEUCHOS_TEST_FOR_EXCEPTION(
480 stepperDIRK == Teuchos::null, std::logic_error,
481 "Error - The implicit component of a general IMEX RK stepper was "
482 "not specified as an DIRK stepper");
483 implicitTableau = stepperDIRK->getTableau();
486 TEUCHOS_TEST_FOR_EXCEPTION(
487 !(explicitTableau != Teuchos::null &&
488 implicitTableau != Teuchos::null),
490 "Error - A parameter list was used to setup a general IMEX RK "
491 "stepper, but did not "
492 "specify both an explicit and an implicit tableau!\n");
494 this->setTableaus(stepperType, explicitTableau, implicitTableau);
495 this->setOrder(pl->get<
int>(
"overall order", 1));
751 this->checkInitialized();
754 using Teuchos::SerialDenseMatrix;
755 using Teuchos::SerialDenseVector;
757 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperIMEX_RK::takeStep()");
759 TEUCHOS_TEST_FOR_EXCEPTION(
760 solutionHistory->getNumStates() < 2, std::logic_error,
761 "Error - StepperIMEX_RK<Scalar>::takeStep(...)\n"
762 <<
"Need at least two SolutionStates for IMEX_RK.\n"
763 <<
" Number of States = " << solutionHistory->getNumStates()
764 <<
"\nTry setting in \"Solution History\" "
765 <<
"\"Storage Type\" = \"Undo\"\n"
766 <<
" or \"Storage Type\" = \"Static\" and "
767 <<
"\"Storage Limit\" = \"2\"\n");
769 RCP<SolutionState<Scalar> > currentState =
770 solutionHistory->getCurrentState();
771 RCP<SolutionState<Scalar> > workingState =
772 solutionHistory->getWorkingState();
773 const Scalar dt = workingState->getTimeStep();
774 const Scalar time = currentState->getTime();
776 const int numStages = explicitTableau_->numStages();
777 const SerialDenseMatrix<int, Scalar>& AHat = explicitTableau_->A();
778 const SerialDenseVector<int, Scalar>& bHat = explicitTableau_->b();
779 const SerialDenseVector<int, Scalar>& cHat = explicitTableau_->c();
780 const SerialDenseMatrix<int, Scalar>& A = implicitTableau_->A();
781 const SerialDenseVector<int, Scalar>& b = implicitTableau_->b();
782 const SerialDenseVector<int, Scalar>& c = implicitTableau_->c();
785 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
787 RCP<StepperIMEX_RK<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
788 this->stepperRKAppAction_->execute(
789 solutionHistory, thisStepper,
793 for (
int i = 0; i < numStages; ++i) {
794 this->setStageNumber(i);
795 Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
796 for (
int j = 0; j < i; ++j) {
797 if (AHat(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
798 Thyra::Vp_StV(xTilde_.ptr(), -dt * AHat(i, j), *(stageF_[j]));
799 if (A(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
800 Thyra::Vp_StV(xTilde_.ptr(), -dt * A(i, j), *(stageG_[j]));
803 this->stepperRKAppAction_->execute(
804 solutionHistory, thisStepper,
807 Scalar ts = time + c(i) * dt;
808 Scalar tHats = time + cHat(i) * dt;
809 if (A(i, i) == Teuchos::ScalarTraits<Scalar>::zero()) {
811 bool isNeeded =
false;
812 for (
int k = i + 1; k < numStages; ++k)
813 if (A(k, i) != 0.0) isNeeded =
true;
814 if (b(i) != 0.0) isNeeded =
true;
815 if (isNeeded ==
false) {
817 assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
820 Thyra::assign(workingState->getX().ptr(), *xTilde_);
821 evalImplicitModelExplicitly(workingState->getX(), ts, dt, i,
827 const Scalar alpha = Scalar(1.0) / (dt * A(i, i));
828 const Scalar beta = Scalar(1.0);
831 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
833 alpha, xTilde_.getConst()));
838 this->stepperRKAppAction_->execute(
839 solutionHistory, thisStepper,
842 const Thyra::SolveStatus<Scalar> sStatus =
843 this->solveImplicitODE(workingState->getX(), stageG_[i], ts, p);
845 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass =
false;
847 this->stepperRKAppAction_->execute(
848 solutionHistory, thisStepper,
852 Thyra::V_StVpStV(stageG_[i].ptr(), -alpha, *workingState->getX(), alpha,
856 this->stepperRKAppAction_->execute(
857 solutionHistory, thisStepper,
859 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
860 this->stepperRKAppAction_->execute(
861 solutionHistory, thisStepper,
866 this->setStageNumber(-1);
869 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
870 for (
int i = 0; i < numStages; ++i) {
871 if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
872 Thyra::Vp_StV((workingState->getX()).ptr(), -dt * bHat(i),
874 if (b(i) != Teuchos::ScalarTraits<Scalar>::zero())
875 Thyra::Vp_StV((workingState->getX()).ptr(), -dt * b(i), *(stageG_[i]));
882 workingState->setOrder(this->getOrder());
883 workingState->computeNorms(currentState);
884 this->stepperRKAppAction_->execute(
885 solutionHistory, thisStepper,