54 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
56 std::string ICConsistency,
57 bool ICConsistencyCheck,
58 bool zeroInitialGuess,
59 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
60 std::string stepperType,
61 Teuchos::RCP<
const RKButcherTableau<Scalar> > explicitTableau,
62 Teuchos::RCP<
const RKButcherTableau<Scalar> > implicitTableau,
65 TEUCHOS_TEST_FOR_EXCEPTION(
66 stepperType !=
"Partitioned IMEX RK 1st order" &&
67 stepperType !=
"Partitioned IMEX RK SSP2" &&
68 stepperType !=
"Partitioned IMEX RK ARS 233" &&
69 stepperType !=
"General Partitioned IMEX RK", std::logic_error,
70 " 'Stepper Type' (='" + stepperType +
"')\n"
71 " does not match one of the types for this Stepper:\n"
72 " 'Partitioned IMEX RK 1st order'\n"
73 " 'Partitioned IMEX RK SSP2'\n"
74 " 'Partitioned IMEX RK ARS 233'\n"
75 " 'General Partitioned IMEX RK'\n");
77 this->setStepperName( stepperType);
78 this->setStepperType( stepperType);
79 this->setUseFSAL( useFSAL);
80 this->setICConsistency( ICConsistency);
81 this->setICConsistencyCheck( ICConsistencyCheck);
82 this->setZeroInitialGuess( zeroInitialGuess);
83 this->setOrder( order);
85 this->setStageNumber(-1);
87 if ( stepperType ==
"General Partitioned IMEX RK" ) {
88 this->setExplicitTableau(explicitTableau);
89 this->setImplicitTableau(implicitTableau);
91 this->setTableaus(stepperType);
93 this->setAppAction(stepperRKAppAction);
94 this->setSolver(solver);
96 if (appModel != Teuchos::null) {
97 this->setModel(appModel);
105 std::string stepperType,
106 Teuchos::RCP<
const RKButcherTableau<Scalar> > explicitTableau,
107 Teuchos::RCP<
const RKButcherTableau<Scalar> > implicitTableau)
109 if (stepperType ==
"") stepperType =
"Partitioned IMEX RK SSP2";
111 if (stepperType ==
"Partitioned IMEX RK 1st order") {
114 typedef Teuchos::ScalarTraits<Scalar> ST;
116 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
117 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
118 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
119 const Scalar one = ST::one();
120 const Scalar zero = ST::zero();
123 A(0,0) = zero; A(0,1) = zero;
124 A(1,0) = one; A(1,1) = zero;
127 b(0) = one; b(1) = zero;
130 c(0) = zero; c(1) = one;
134 auto expTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
135 "Explicit Tableau - Partitioned IMEX RK 1st order",
136 A,b,c,order,order,order));
137 expTableau->setTVD(
true);
138 expTableau->setTVDCoeff(2.0);
140 this->setExplicitTableau(expTableau);
144 typedef Teuchos::ScalarTraits<Scalar> ST;
146 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
147 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
148 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
149 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
150 const Scalar one = ST::one();
151 const Scalar zero = ST::zero();
154 A(0,0) = zero; A(0,1) = zero;
155 A(1,0) = zero; A(1,1) = one;
158 b(0) = zero; b(1) = one;
161 c(0) = zero; c(1) = one;
165 auto impTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
166 "Implicit Tableau - Partitioned IMEX RK 1st order",
167 A,b,c,order,order,order));
168 impTableau->setTVD(
true);
169 impTableau->setTVDCoeff(sspcoef);
171 this->setImplicitTableau(impTableau);
173 this->setStepperName(
"Partitioned IMEX RK 1st order");
174 this->setStepperType(
"Partitioned IMEX RK 1st order");
177 }
else if (stepperType ==
"Partitioned IMEX RK SSP2") {
179 auto stepperERK = Teuchos::rcp(
new StepperERK_Trapezoidal<Scalar>());
180 this->setExplicitTableau(stepperERK->getTableau());
183 auto stepperSDIRK = Teuchos::rcp(
new StepperSDIRK_2Stage3rdOrder<Scalar>());
184 stepperSDIRK->setGammaType(
"2nd Order L-stable");
185 this->setImplicitTableau(stepperSDIRK->getTableau());
187 this->setStepperName(
"Partitioned IMEX RK SSP2");
188 this->setStepperType(
"Partitioned IMEX RK SSP2");
191 }
else if (stepperType ==
"Partitioned IMEX RK ARS 233") {
192 typedef Teuchos::ScalarTraits<Scalar> ST;
194 Teuchos::SerialDenseMatrix<int,Scalar> A(NumStages,NumStages);
195 Teuchos::SerialDenseVector<int,Scalar> b(NumStages);
196 Teuchos::SerialDenseVector<int,Scalar> c(NumStages);
197 const Scalar one = ST::one();
198 const Scalar zero = ST::zero();
199 const Scalar onehalf = ST::one()/(2*ST::one());
200 const Scalar gamma = (3*one+ST::squareroot(3*one))/(6*one);
204 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
205 A(1,0) = gamma; A(1,1) = zero; A(1,2) = zero;
206 A(2,0) = (gamma-1.0); A(2,1) = (2.0-2.0*gamma); A(2,2) = zero;
209 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
212 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
216 auto expTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
217 "Explicit Tableau - Partitioned IMEX RK ARS 233",
218 A,b,c,order,order,order));
220 this->setExplicitTableau(expTableau);
225 A(0,0) = zero; A(0,1) = zero; A(0,2) = zero;
226 A(1,0) = zero; A(1,1) = gamma; A(1,2) = zero;
227 A(2,0) = zero; A(2,1) = (1.0-2.0*gamma); A(2,2) = gamma;
230 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
233 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
237 auto impTableau = Teuchos::rcp(
new RKButcherTableau<Scalar>(
238 "Implicit Tableau - Partitioned IMEX RK ARS 233",
239 A,b,c,order,order,order));
241 this->setImplicitTableau(impTableau);
243 this->setStepperName(
"Partitioned IMEX RK ARS 233");
244 this->setStepperType(
"Partitioned IMEX RK ARS 233");
247 }
else if (stepperType ==
"General Partitioned IMEX RK") {
249 if ( explicitTableau == Teuchos::null ) {
251 auto stepperERK = Teuchos::rcp(
new StepperERK_Trapezoidal<Scalar>());
252 this->setExplicitTableau(stepperERK->getTableau());
254 this->setExplicitTableau(explicitTableau);
257 if ( implicitTableau == Teuchos::null ) {
259 auto stepperSDIRK = Teuchos::rcp(
new StepperSDIRK_2Stage3rdOrder<Scalar>());
260 stepperSDIRK->setGammaType(
"2nd Order L-stable");
261 this->setImplicitTableau(stepperSDIRK->getTableau());
263 this->setImplicitTableau(implicitTableau);
266 this->setStepperName(
"General Partitioned IMEX RK");
267 this->setStepperType(
"General Partitioned IMEX RK");
271 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
272 "Error - Not a valid StepperIMEX_RK_Partition type! Stepper Type = "
273 << stepperType <<
"\n"
274 <<
" Current valid types are: " <<
"\n"
275 <<
" 'Partitioned IMEX RK 1st order'" <<
"\n"
276 <<
" 'Partitioned IMEX RK SSP2'" <<
"\n"
277 <<
" 'Partitioned IMEX RK ARS 233'" <<
"\n"
278 <<
" 'General Partitioned IMEX RK'" <<
"\n");
281 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_==Teuchos::null,
283 "Error - StepperIMEX_RK_Partition - Explicit tableau is null!");
284 TEUCHOS_TEST_FOR_EXCEPTION(implicitTableau_==Teuchos::null,
286 "Error - StepperIMEX_RK_Partition - Implicit tableau is null!");
287 TEUCHOS_TEST_FOR_EXCEPTION(
288 explicitTableau_->numStages()!=implicitTableau_->numStages(),
290 "Error - StepperIMEX_RK_Partition - Number of stages do not match!\n"
291 <<
" Explicit tableau = " << explicitTableau_->description() <<
"\n"
292 <<
" number of stages = " << explicitTableau_->numStages() <<
"\n"
293 <<
" Implicit tableau = " << implicitTableau_->description() <<
"\n"
294 <<
" number of stages = " << implicitTableau_->numStages() <<
"\n");
296 this->isInitialized_ =
false;
304 Teuchos::RCP<Teuchos::ParameterList> pl,
305 std::string stepperType)
308 if (stepperType ==
"") {
309 if (pl == Teuchos::null)
310 stepperType =
"Partitioned IMEX RK SSP2";
312 stepperType = pl->get<std::string>(
"Stepper Type",
"Partitioned IMEX RK SSP2");
315 if (stepperType !=
"General Partitioned IMEX RK") {
316 this->setTableaus(stepperType);
318 if (pl != Teuchos::null) {
319 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau;
320 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau;
321 if (pl->isSublist(
"IMEX-RK Explicit Stepper")) {
322 RCP<Teuchos::ParameterList> explicitPL = Teuchos::rcp(
323 new Teuchos::ParameterList(pl->sublist(
"IMEX-RK Explicit Stepper")));
324 auto sf = Teuchos::rcp(
new StepperFactory<Scalar>());
325 auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
326 auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> > (
328 TEUCHOS_TEST_FOR_EXCEPTION(stepperERK == Teuchos::null, std::logic_error,
329 "Error - The explicit component of a general partitioned IMEX RK stepper was not specified as an ExplicitRK stepper");
330 explicitTableau = stepperERK->getTableau();
333 if (pl->isSublist(
"IMEX-RK Implicit Stepper")) {
334 RCP<Teuchos::ParameterList> implicitPL = Teuchos::rcp(
335 new Teuchos::ParameterList(pl->sublist(
"IMEX-RK Implicit Stepper")));
336 auto sf = Teuchos::rcp(
new StepperFactory<Scalar>());
337 auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
338 auto stepperDIRK = Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> > (
340 TEUCHOS_TEST_FOR_EXCEPTION(stepperDIRK == Teuchos::null, std::logic_error,
341 "Error - The implicit component of a general partitioned IMEX RK stepper was not specified as an DIRK stepper");
342 implicitTableau = stepperDIRK->getTableau();
345 TEUCHOS_TEST_FOR_EXCEPTION(
346 !(explicitTableau!=Teuchos::null && implicitTableau!=Teuchos::null), std::logic_error,
347 "Error - A parameter list was used to setup a general partitioned IMEX RK stepper, but did not "
348 "specify both an explicit and an implicit tableau!\n");
350 this->setTableaus(stepperType, explicitTableau, implicitTableau);
352 this->setOrder(pl->get<
int>(
"overall order", 1));
606 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
608 this->checkInitialized();
611 using Teuchos::SerialDenseMatrix;
612 using Teuchos::SerialDenseVector;
614 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperIMEX_RK_Partition::takeStep()");
616 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
618 "Error - StepperIMEX_RK_Partition<Scalar>::takeStep(...)\n"
619 "Need at least two SolutionStates for IMEX_RK_Partition.\n"
620 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
621 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
622 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
624 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
625 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
626 const Scalar dt = workingState->getTimeStep();
627 const Scalar time = currentState->getTime();
629 const int numStages = explicitTableau_->numStages();
630 const SerialDenseMatrix<int,Scalar> & AHat = explicitTableau_->A();
631 const SerialDenseVector<int,Scalar> & bHat = explicitTableau_->b();
632 const SerialDenseVector<int,Scalar> & cHat = explicitTableau_->c();
633 const SerialDenseMatrix<int,Scalar> & A = implicitTableau_->A();
634 const SerialDenseVector<int,Scalar> & b = implicitTableau_->b();
635 const SerialDenseVector<int,Scalar> & c = implicitTableau_->c();
637 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
638 wrapperModelPairIMEX =
639 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
640 (this->wrapperModel_);
643 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
644 RCP<Thyra::VectorBase<Scalar> > stageY =
645 wrapperModelPairIMEX->getExplicitOnlyVector(workingState->getX());
646 RCP<Thyra::VectorBase<Scalar> > stageX =
647 wrapperModelPairIMEX->getIMEXVector(workingState->getX());
649 RCP<StepperIMEX_RK_Partition<Scalar> > thisStepper = Teuchos::rcpFromRef(*
this);
650 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
654 for (
int i = 0; i < numStages; ++i) {
655 this->setStageNumber(i);
657 Thyra::assign(stageY.ptr(),
658 *(wrapperModelPairIMEX->getExplicitOnlyVector(currentState->getX())));
659 Thyra::assign(xTilde_.ptr(),
660 *(wrapperModelPairIMEX->getIMEXVector(currentState->getX())));
661 for (
int j = 0; j < i; ++j) {
662 if (AHat(i,j) != Teuchos::ScalarTraits<Scalar>::zero()) {
663 RCP<Thyra::VectorBase<Scalar> > stageFy =
664 wrapperModelPairIMEX->getExplicitOnlyVector(stageF_[j]);
665 RCP<Thyra::VectorBase<Scalar> > stageFx =
666 wrapperModelPairIMEX->getIMEXVector(stageF_[j]);
667 Thyra::Vp_StV(stageY.ptr(), -dt*AHat(i,j), *stageFy);
668 Thyra::Vp_StV(xTilde_.ptr(), -dt*AHat(i,j), *stageFx);
670 if (A (i,j) != Teuchos::ScalarTraits<Scalar>::zero())
671 Thyra::Vp_StV(xTilde_.ptr(), -dt*A (i,j), *(stageGx_[j]));
674 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
677 Scalar ts = time + c(i)*dt;
678 Scalar tHats = time + cHat(i)*dt;
679 if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
681 bool isNeeded =
false;
682 for (
int k=i+1; k<numStages; ++k)
if (A(k,i) != 0.0) isNeeded =
true;
683 if (b(i) != 0.0) isNeeded =
true;
684 if (isNeeded ==
false) {
686 assign(stageGx_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
688 Thyra::assign(stageX.ptr(), *xTilde_);
689 evalImplicitModelExplicitly(stageX, stageY, ts, dt, i, stageGx_[i]);
693 const Scalar alpha = Scalar(1.0)/(dt*A(i,i));
694 const Scalar beta = Scalar(1.0);
697 RCP<TimeDerivative<Scalar> > timeDer =
698 Teuchos::rcp(
new StepperIMEX_RKPartTimeDerivative<Scalar>(
699 alpha, xTilde_.getConst()));
701 auto p = Teuchos::rcp(
new ImplicitODEParameters<Scalar>(
704 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
707 wrapperModelPairIMEX->setUseImplicitModel(
true);
708 this->solver_->setModel(wrapperModelPairIMEX);
710 Thyra::SolveStatus<Scalar> sStatus;
711 if (wrapperModelPairIMEX->getParameterIndex() >= 0)
712 sStatus = this->solveImplicitODE(stageX, stageGx_[i], ts, p,
713 stageY, wrapperModelPairIMEX->getParameterIndex());
715 sStatus = this->solveImplicitODE(stageX, stageGx_[i], ts, p);
717 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass =
false;
719 wrapperModelPairIMEX->setUseImplicitModel(
false);
721 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
725 Thyra::V_StVpStV(stageGx_[i].ptr(), -alpha, *stageX, alpha, *xTilde_);
728 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
730 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
731 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
737 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
738 RCP<Thyra::VectorBase<Scalar> > Z = workingState->getX();
739 RCP<Thyra::VectorBase<Scalar> > X = wrapperModelPairIMEX->getIMEXVector(Z);
740 for (
int i=0; i < numStages; ++i) {
741 if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
742 Thyra::Vp_StV(Z.ptr(), -dt*bHat(i), *(stageF_[i]));
743 if (b (i) != Teuchos::ScalarTraits<Scalar>::zero())
744 Thyra::Vp_StV(X.ptr(), -dt*b (i), *(stageGx_[i]));
747 if (pass ==
true) workingState->setSolutionStatus(
Status::PASSED);
749 workingState->setOrder(this->getOrder());
750 workingState->computeNorms(currentState);
751 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
755 this->setStageNumber(-1);
864 auto pl = this->getValidParametersBasicImplicit();
865 pl->template set<int>(
"overall order", this->getOrder());
867 auto explicitStepper = Teuchos::rcp(
new StepperERK_General<Scalar>());
868 explicitStepper->setTableau(
869 explicitTableau_->A(), explicitTableau_->b(), explicitTableau_->c(),
870 explicitTableau_->order(), explicitTableau_->orderMin(),
871 explicitTableau_->orderMax(), explicitTableau_->bstar() );
872 pl->set(
"IMEX-RK Explicit Stepper", *explicitStepper->getValidParameters());
874 auto implicitStepper = Teuchos::rcp(
new StepperERK_General<Scalar>());
875 implicitStepper->setTableau(
876 implicitTableau_->A(), implicitTableau_->b(), implicitTableau_->c(),
877 implicitTableau_->order(), implicitTableau_->orderMin(),
878 implicitTableau_->orderMax(), implicitTableau_->bstar() );
879 pl->set(
"IMEX-RK Implicit Stepper", *implicitStepper->getValidParameters());