133 Teuchos::RCP<
const RKButcherTableau<Scalar> > explicitTableau,
134 Teuchos::RCP<
const RKButcherTableau<Scalar> > implicitTableau)
136 if (stepperType ==
"") stepperType =
"IMEX RK SSP2";
138 if (stepperType ==
"IMEX RK 1st order" ) {
141 typedef Teuchos::ScalarTraits<Scalar> ST;
143 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
144 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
145 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
146 const Scalar one = ST::one();
147 const Scalar zero = ST::zero();
150 A(0,0) = zero; A(0,1) = zero;
151 A(1,0) = one; A(1,1) = zero;
154 b(0) = one; b(1) = zero;
157 c(0) = zero; c(1) = one;
161 auto expTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
162 "Explicit Tableau - IMEX RK 1st order",
163 A,b,c,order,order,order));
164 expTableau->setTVD(
true);
165 expTableau->setTVDCoeff(2.0);
167 this->setExplicitTableau(expTableau);
171 typedef Teuchos::ScalarTraits<Scalar> ST;
173 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
174 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
175 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
176 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
177 const Scalar one = ST::one();
178 const Scalar zero = ST::zero();
181 A(0,0) = zero; A(0,1) = zero;
182 A(1,0) = zero; A(1,1) = one;
185 b(0) = zero; b(1) = one;
188 c(0) = zero; c(1) = one;
192 auto impTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
193 "Implicit Tableau - IMEX RK 1st order",
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" )
208 typedef Teuchos::ScalarTraits<Scalar> ST;
209 const int NumStages = 1;
211 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
212 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
213 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
214 const Scalar one = ST::one();
215 const Scalar zero = ST::zero();
226 auto expTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
227 "Explicit Tableau - SSP1_111",A,b,c,order,order,order));
228 expTableau->setTVD(
true);
229 expTableau->setTVDCoeff(1.0);
231 this->setExplicitTableau(expTableau);
235 typedef Teuchos::ScalarTraits<Scalar> ST;
236 const int NumStages = 1;
238 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
239 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
240 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
241 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
242 const Scalar one = ST::one();
253 auto impTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
254 "Implicit Tableau - SSP1_111",
255 A,b,c,order,order,order));
256 impTableau->setTVD(
true);
257 impTableau->setTVDCoeff(sspcoef);
259 this->setImplicitTableau(impTableau);
261 this->setStepperName(
"SSP1_111");
262 this->setStepperType(
"SSP1_111");
265 }
else if (stepperType ==
"IMEX RK SSP2" || stepperType ==
"SSP2_222_L" ) {
267 auto stepperERK = Teuchos::rcp(
new StepperERK_Trapezoidal<Scalar>());
268 this->setExplicitTableau(stepperERK->getTableau());
271 auto stepperSDIRK = Teuchos::rcp(
new StepperSDIRK_2Stage3rdOrder<Scalar>());
272 stepperSDIRK->setGammaType(
"2nd Order L-stable");
273 this->setImplicitTableau(stepperSDIRK->getTableau());
275 this->setStepperName(
"IMEX RK SSP2");
276 this->setStepperType(
"IMEX RK SSP2");
278 }
else if (stepperType ==
"SSP2_222" || stepperType ==
"SSP2_222_A" ) {
280 auto stepperERK = Teuchos::rcp(
new StepperERK_Trapezoidal<Scalar>());
281 this->setExplicitTableau(stepperERK->getTableau());
284 auto stepperSDIRK = Teuchos::rcp(
new StepperSDIRK_2Stage3rdOrder<Scalar>());
285 stepperSDIRK->setGammaType(
"gamma");
286 stepperSDIRK->setGamma( 0.5 );
287 this->setImplicitTableau(stepperSDIRK->getTableau());
289 this->setStepperName(
"SSP2_222");
290 this->setStepperType(
"SSP2_222");
292 }
else if (stepperType ==
"IMEX RK SSP3" || stepperType ==
"SSP3_332" ) {
294 auto stepperERK = Teuchos::rcp(
new StepperERK_3Stage3rdOrderTVD<Scalar>());
295 this->setExplicitTableau(stepperERK->getTableau());
298 auto stepperSDIRK = Teuchos::rcp(
new StepperSDIRK_3Stage2ndOrder<Scalar>());
299 this->setImplicitTableau(stepperSDIRK->getTableau());
301 this->setStepperName(
"IMEX RK SSP3");
302 this->setStepperType(
"IMEX RK SSP3");
304 }
else if (stepperType ==
"IMEX RK ARS 233" || stepperType ==
"ARS 233" ) {
305 typedef Teuchos::ScalarTraits<Scalar> ST;
307 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
308 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
309 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
310 const Scalar one = ST::one();
311 const Scalar zero = ST::zero();
312 const Scalar onehalf = ST::one()/(2*ST::one());
313 const Scalar gamma = (3*one+ST::squareroot(3*one))/(6*one);
317 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
318 A(1,0) = gamma; A(1,1) = zero; A(1,2) = zero;
319 A(2,0) = (gamma-1.0); A(2,1) = (2.0-2.0*gamma); A(2,2) = zero;
322 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
325 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
329 auto expTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
330 "Partition IMEX-RK Explicit Stepper",A,b,c,order,order,order));
332 this->setExplicitTableau(expTableau);
337 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
338 A(1,0) = zero; A(1,1) = gamma; A(1,2) = zero;
339 A(2,0) = zero; A(2,1) = (1.0-2.0*gamma); A(2,2) = gamma;
342 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
345 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
349 auto impTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
350 "Partition IMEX-RK Implicit Stepper",A,b,c,order,order,order));
352 this->setImplicitTableau(impTableau);
354 this->setStepperName(
"IMEX RK ARS 233");
355 this->setStepperType(
"IMEX RK ARS 233");
358 }
else if (stepperType ==
"General IMEX RK") {
360 if (explicitTableau == Teuchos::null) {
362 auto stepperERK = Teuchos::rcp(
new StepperERK_Trapezoidal<Scalar>());
363 this->setExplicitTableau(stepperERK->getTableau());
365 this->setExplicitTableau(explicitTableau);
368 if (explicitTableau == Teuchos::null) {
370 auto stepperSDIRK = Teuchos::rcp(
new StepperSDIRK_2Stage3rdOrder<Scalar>());
371 stepperSDIRK->setGammaType(
"2nd Order L-stable");
372 this->setImplicitTableau(stepperSDIRK->getTableau());
374 this->setImplicitTableau(implicitTableau);
377 this->setStepperName(
"General IMEX RK");
378 this->setStepperType(
"General IMEX RK");
382 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
383 "Error - Not a valid StepperIMEX_RK type! Stepper Type = "
384 << stepperType <<
"\n"
385 <<
" Current valid types are: \n"
386 <<
" 'IMEX RK 1st order\n"
388 <<
" 'IMEX RK SSP2' ('SSP2_222_L')\n"
389 <<
" 'SSP2_222' ('SSP2_222_A')\n"
390 <<
" 'IMEX RK SSP3' ('SSP3_332')\n"
391 <<
" 'IMEX RK ARS 233' ('ARS 233')\n"
392 <<
" 'General IMEX RK'\n");
395 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_==Teuchos::null,
397 "Error - StepperIMEX_RK - Explicit tableau is null!");
398 TEUCHOS_TEST_FOR_EXCEPTION(implicitTableau_==Teuchos::null,
400 "Error - StepperIMEX_RK - Implicit tableau is null!");
401 TEUCHOS_TEST_FOR_EXCEPTION(
402 explicitTableau_->numStages()!=implicitTableau_->numStages(),
404 "Error - StepperIMEX_RK - Number of stages do not match!\n"
405 <<
" Explicit tableau = " << explicitTableau_->description() <<
"\n"
406 <<
" number of stages = " << explicitTableau_->numStages() <<
"\n"
407 <<
" Implicit tableau = " << implicitTableau_->description() <<
"\n"
408 <<
" number of stages = " << implicitTableau_->numStages() <<
"\n");
410 this->isInitialized_ =
false;
417setTableaus(Teuchos::RCP<Teuchos::ParameterList> pl,
418 std::string stepperType)
421 if (stepperType ==
"") {
422 if (pl == Teuchos::null)
423 stepperType =
"IMEX RK SSP2";
425 stepperType = pl->get<std::string>(
"Stepper Type",
"IMEX RK SSP2");
428 if (stepperType !=
"General IMEX RK") {
429 this->setTableaus(stepperType);
431 if (pl != Teuchos::null) {
432 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau;
433 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau;
434 if (pl->isSublist(
"IMEX-RK Explicit Stepper")) {
435 RCP<Teuchos::ParameterList> explicitPL = Teuchos::rcp(
436 new Teuchos::ParameterList(pl->sublist(
"IMEX-RK Explicit Stepper")));
437 auto sf = Teuchos::rcp(
new StepperFactory<Scalar>());
438 auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
439 auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> > (
441 TEUCHOS_TEST_FOR_EXCEPTION(stepperERK == Teuchos::null, std::logic_error,
442 "Error - The explicit component of a general IMEX RK stepper was not specified as an ExplicitRK stepper");
443 explicitTableau = stepperERK->getTableau();
446 if (pl->isSublist(
"IMEX-RK Implicit Stepper")) {
447 RCP<Teuchos::ParameterList> implicitPL = Teuchos::rcp(
448 new Teuchos::ParameterList(pl->sublist(
"IMEX-RK Implicit Stepper")));
449 auto sf = Teuchos::rcp(
new StepperFactory<Scalar>());
450 auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
451 auto stepperDIRK = Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> > (
453 TEUCHOS_TEST_FOR_EXCEPTION(stepperDIRK == Teuchos::null, std::logic_error,
454 "Error - The implicit component of a general IMEX RK stepper was not specified as an DIRK stepper");
455 implicitTableau = stepperDIRK->getTableau();
458 TEUCHOS_TEST_FOR_EXCEPTION(
459 !(explicitTableau!=Teuchos::null && implicitTableau!=Teuchos::null), std::logic_error,
460 "Error - A parameter list was used to setup a general IMEX RK stepper, but did not "
461 "specify both an explicit and an implicit tableau!\n");
463 this->setTableaus(stepperType, explicitTableau, implicitTableau);
464 this->setOrder(pl->get<
int>(
"overall order", 1));
717 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
719 this->checkInitialized();
722 using Teuchos::SerialDenseMatrix;
723 using Teuchos::SerialDenseVector;
725 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperIMEX_RK::takeStep()");
727 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
729 "Error - StepperIMEX_RK<Scalar>::takeStep(...)\n"
730 "Need at least two SolutionStates for IMEX_RK.\n"
731 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
732 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
733 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
735 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
736 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
737 const Scalar dt = workingState->getTimeStep();
738 const Scalar time = currentState->getTime();
740 const int numStages = explicitTableau_->numStages();
741 const SerialDenseMatrix<int,Scalar> & AHat = explicitTableau_->A();
742 const SerialDenseVector<int,Scalar> & bHat = explicitTableau_->b();
743 const SerialDenseVector<int,Scalar> & cHat = explicitTableau_->c();
744 const SerialDenseMatrix<int,Scalar> & A = implicitTableau_->A();
745 const SerialDenseVector<int,Scalar> & b = implicitTableau_->b();
746 const SerialDenseVector<int,Scalar> & c = implicitTableau_->c();
749 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
751 RCP<StepperIMEX_RK<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
752 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
756 for (
int i = 0; i < numStages; ++i) {
757 this->setStageNumber(i);
758 Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
759 for (
int j = 0; j < i; ++j) {
760 if (AHat(i,j) != Teuchos::ScalarTraits<Scalar>::zero())
761 Thyra::Vp_StV(xTilde_.ptr(), -dt*AHat(i,j), *(stageF_[j]));
762 if (A (i,j) != Teuchos::ScalarTraits<Scalar>::zero())
763 Thyra::Vp_StV(xTilde_.ptr(), -dt*A (i,j), *(stageG_[j]));
766 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
769 Scalar ts = time + c(i)*dt;
770 Scalar tHats = time + cHat(i)*dt;
771 if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
773 bool isNeeded =
false;
774 for (
int k=i+1; k<numStages; ++k)
if (A(k,i) != 0.0) isNeeded =
true;
775 if (b(i) != 0.0) isNeeded =
true;
776 if (isNeeded ==
false) {
778 assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
780 Thyra::assign(workingState->getX().ptr(), *xTilde_);
781 evalImplicitModelExplicitly(workingState->getX(), ts, dt, i, stageG_[i]);
785 const Scalar alpha = Scalar(1.0)/(dt*A(i,i));
786 const Scalar beta = Scalar(1.0);
789 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
790 Teuchos::rcp(
new StepperIMEX_RKTimeDerivative<Scalar>(
791 alpha, xTilde_.getConst()));
793 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
796 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
799 const Thyra::SolveStatus<Scalar> sStatus =
800 this->solveImplicitODE(workingState->getX(), stageG_[i], ts, p);
802 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass =
false;
804 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
808 Thyra::V_StVpStV(stageG_[i].ptr(), -alpha, *workingState->getX(), alpha, *xTilde_);
811 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
813 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
814 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
819 this->setStageNumber(-1);
822 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
823 for (
int i=0; i < numStages; ++i) {
824 if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
825 Thyra::Vp_StV((workingState->getX()).ptr(), -dt*bHat(i), *(stageF_[i]));
826 if (b (i) != Teuchos::ScalarTraits<Scalar>::zero())
827 Thyra::Vp_StV((workingState->getX()).ptr(), -dt*b (i), *(stageG_[i]));
830 if (pass ==
true) workingState->setSolutionStatus(
Status::PASSED);
832 workingState->setOrder(this->getOrder());
833 workingState->computeNorms(currentState);
834 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,