9#ifndef Tempus_StepperIMEX_RK_impl_hpp
10#define Tempus_StepperIMEX_RK_impl_hpp
12#include "Thyra_VectorStdOps.hpp"
14#include "Tempus_StepperFactory.hpp"
24 TEUCHOS_TEST_FOR_EXCEPTION(
25 stepperType !=
"IMEX RK 1st order" &&
26 stepperType !=
"SSP1_111" &&
27 stepperType !=
"IMEX RK SSP2" &&
28 stepperType !=
"IMEX RK SSP3" &&
29 stepperType !=
"SSP3_332" &&
30 stepperType !=
"SSP2_222" &&
31 stepperType !=
"SSP2_222_L" &&
32 stepperType !=
"SSP2_222_A" &&
33 stepperType !=
"IMEX RK ARS 233" &&
34 stepperType !=
"ARS 233" &&
35 stepperType !=
"General IMEX RK", std::logic_error,
36 " 'Stepper Type' (='" + stepperType +
"')\n"
37 " does not match one of the types for this Stepper:\n"
38 " 'IMEX RK 1st order'\n"
46 " 'IMEX RK ARS 233'\n"
48 " 'General IMEX RK'\n");
50 this->setStepperName( stepperType);
51 this->setStepperType( stepperType);
52 this->setUseFSAL(
false);
53 this->setICConsistency(
"None");
54 this->setICConsistencyCheck(
false);
55 this->setZeroInitialGuess(
false);
57 this->setStageNumber(-1);
59 this->setTableaus(stepperType);
60 this->setAppAction(Teuchos::null);
61 this->setDefaultSolver();
68 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
70 std::string ICConsistency,
71 bool ICConsistencyCheck,
72 bool zeroInitialGuess,
74 std::string stepperType,
79 TEUCHOS_TEST_FOR_EXCEPTION(
80 stepperType !=
"IMEX RK 1st order" &&
81 stepperType !=
"SSP1_111" &&
82 stepperType !=
"IMEX RK SSP2" &&
83 stepperType !=
"IMEX RK SSP3" &&
84 stepperType !=
"SSP3_332" &&
85 stepperType !=
"SSP2_222" &&
86 stepperType !=
"SSP2_222_L" &&
87 stepperType !=
"SSP2_222_A" &&
88 stepperType !=
"IMEX RK ARS 233" &&
89 stepperType !=
"ARS 233" &&
90 stepperType !=
"General IMEX RK", std::logic_error,
91 " 'Stepper Type' (='" + stepperType +
"')\n"
92 " does not match one of the types for this Stepper:\n"
93 " 'IMEX RK 1st order'\n"
101 " 'IMEX RK ARS 233'\n"
103 " 'General IMEX RK'\n");
105 this->setStepperName( stepperType);
106 this->setStepperType( stepperType);
107 this->setUseFSAL( useFSAL);
108 this->setICConsistency( ICConsistency);
109 this->setICConsistencyCheck( ICConsistencyCheck);
110 this->setZeroInitialGuess( zeroInitialGuess);
112 this->setStageNumber(-1);
114 if ( stepperType ==
"General IMEX RK" ) {
115 this->setExplicitTableau(explicitTableau);
116 this->setImplicitTableau(implicitTableau);
118 this->setTableaus(stepperType);
120 this->setOrder(order);
121 this->setAppAction(stepperRKAppAction);
122 this->setSolver(solver);
124 if (appModel != Teuchos::null) {
125 this->setModel(appModel);
131template<
class Scalar>
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;
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;
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();
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();
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" ) {
268 this->setExplicitTableau(stepperERK->getTableau());
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" ) {
281 this->setExplicitTableau(stepperERK->getTableau());
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" ) {
295 this->setExplicitTableau(stepperERK->getTableau());
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;
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;
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) {
363 this->setExplicitTableau(stepperERK->getTableau());
365 this->setExplicitTableau(explicitTableau);
368 if (explicitTableau == Teuchos::null) {
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;
414template<
class Scalar>
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")));
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")));
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));
470template<
class Scalar>
474 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau->isImplicit() ==
true,
476 "Error - Received an implicit Tableau for setExplicitTableau()!\n" <<
477 " Tableau = " << explicitTableau->description() <<
"\n");
478 explicitTableau_ = explicitTableau;
480 this->isInitialized_ =
false;
484template<
class Scalar>
488 TEUCHOS_TEST_FOR_EXCEPTION( implicitTableau->isDIRK() !=
true,
490 "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n" <<
491 " Tableau = " << implicitTableau->description() <<
"\n");
492 implicitTableau_ = implicitTableau;
494 this->isInitialized_ =
false;
497template<
class Scalar>
502 using Teuchos::rcp_const_cast;
503 using Teuchos::rcp_dynamic_cast;
504 RCP<Thyra::ModelEvaluator<Scalar> > ncModel =
505 rcp_const_cast<Thyra::ModelEvaluator<Scalar> > (appModel);
506 RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > modelPairIMEX =
507 rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > (ncModel);
508 TEUCHOS_TEST_FOR_EXCEPTION( modelPairIMEX == Teuchos::null, std::logic_error,
509 "Error - StepperIMEX_RK::setModel() was given a ModelEvaluator that\n"
510 " could not be cast to a WrapperModelEvaluatorPairIMEX_Basic!\n"
511 " From: " << appModel <<
"\n"
512 " To : " << modelPairIMEX <<
"\n"
513 " Likely have given the wrong ModelEvaluator to this Stepper.\n");
515 setModelPair(modelPairIMEX);
517 TEUCHOS_TEST_FOR_EXCEPTION(this->solver_ == Teuchos::null, std::logic_error,
518 "Error - Solver is not set!\n");
519 this->solver_->setModel(modelPairIMEX);
521 this->isInitialized_ =
false;
530template<
class Scalar>
535 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
536 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
537 this->wrapperModel_);
540 wrapperModelPairIMEX = modelPairIMEX;
541 wrapperModelPairIMEX->initialize();
542 int expXDim = wrapperModelPairIMEX->getExplicitModel()->get_x_space()->dim();
543 int impXDim = wrapperModelPairIMEX->getImplicitModel()->get_x_space()->dim();
544 TEUCHOS_TEST_FOR_EXCEPTION( expXDim != impXDim, std::logic_error,
546 " Explicit and Implicit x vectors are incompatible!\n"
547 " Explicit vector dim = " << expXDim <<
"\n"
548 " Implicit vector dim = " << impXDim <<
"\n");
550 this->wrapperModel_ = wrapperModelPairIMEX;
552 this->isInitialized_ =
false;
560template<
class Scalar>
567 this->wrapperModel_ = Teuchos::rcp(
569 explicitModel, implicitModel));
571 this->isInitialized_ =
false;
575template<
class Scalar>
578 TEUCHOS_TEST_FOR_EXCEPTION( this->wrapperModel_ == Teuchos::null,
580 "Error - Need to set the model, setModel(), before calling "
581 "StepperIMEX_RK::initialize()\n");
584 const int numStages = explicitTableau_->numStages();
585 stageF_.resize(numStages);
586 stageG_.resize(numStages);
587 for(
int i=0; i < numStages; i++) {
588 stageF_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
589 stageG_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
590 assign(stageF_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
591 assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
594 xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
595 assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
601template<
class Scalar>
607 int numStates = solutionHistory->getNumStates();
609 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
610 "Error - setInitialConditions() needs at least one SolutionState\n"
611 " to set the initial condition. Number of States = " << numStates);
614 RCP<Teuchos::FancyOStream> out = this->getOStream();
615 Teuchos::OSTab ostab(out,1,
"StepperIMEX_RK::setInitialConditions()");
616 *out <<
"Warning -- SolutionHistory has more than one state!\n"
617 <<
"Setting the initial conditions on the currentState.\n"<<std::endl;
620 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
621 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
624 auto inArgs = this->wrapperModel_->getNominalValues();
625 if (x == Teuchos::null) {
626 TEUCHOS_TEST_FOR_EXCEPTION( (x == Teuchos::null) &&
627 (inArgs.get_x() == Teuchos::null), std::logic_error,
628 "Error - setInitialConditions() needs the ICs from the SolutionHistory\n"
629 " or getNominalValues()!\n");
631 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
632 initialState->setX(x);
636 std::string icConsistency = this->getICConsistency();
637 TEUCHOS_TEST_FOR_EXCEPTION(icConsistency !=
"None", std::logic_error,
638 "Error - setInitialConditions() requested a consistency of '"
639 << icConsistency <<
"'.\n"
640 " But only 'None' is available for IMEX-RK!\n");
642 TEUCHOS_TEST_FOR_EXCEPTION( this->getUseFSAL(), std::logic_error,
643 "Error - The First-Same-As-Last (FSAL) principle is not "
644 <<
"available for IMEX-RK. Set useFSAL=false.\n");
648template <
typename Scalar>
651 Scalar time, Scalar stepSize, Scalar stageNumber,
654 typedef Thyra::ModelEvaluatorBase MEB;
655 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
656 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
657 this->wrapperModel_);
658 MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->getInArgs();
660 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
661 if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
662 if (inArgs.supports(MEB::IN_ARG_stage_number))
663 inArgs.set_stage_number(stageNumber);
670 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
672 MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->getOutArgs();
675 wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs,outArgs);
676 Thyra::Vt_S(G.ptr(), -1.0);
680template <
typename Scalar>
683 Scalar time, Scalar stepSize, Scalar stageNumber,
686 typedef Thyra::ModelEvaluatorBase MEB;
688 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
689 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
690 this->wrapperModel_);
691 MEB::InArgs<Scalar> inArgs =
692 wrapperModelPairIMEX->getExplicitModel()->createInArgs();
694 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
695 if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
696 if (inArgs.supports(MEB::IN_ARG_stage_number))
697 inArgs.set_stage_number(stageNumber);
704 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
706 MEB::OutArgs<Scalar> outArgs =
707 wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
710 wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
711 Thyra::Vt_S(F.ptr(), -1.0);
715template<
class Scalar>
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 =
791 alpha, xTilde_.getConst()));
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,
846template<
class Scalar>
847Teuchos::RCP<Tempus::StepperState<Scalar> >
851 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
857template<
class Scalar>
859 Teuchos::FancyOStream &out,
860 const Teuchos::EVerbosityLevel verbLevel)
const
862 out.setOutputToRootOnly(0);
868 out <<
"--- StepperIMEX_RK_Partition ---\n";
869 out <<
" explicitTableau_ = " << explicitTableau_ << std::endl;
870 if (verbLevel == Teuchos::VERB_HIGH)
871 explicitTableau_->describe(out, verbLevel);
872 out <<
" implicitTableau_ = " << implicitTableau_ << std::endl;
873 if (verbLevel == Teuchos::VERB_HIGH)
874 implicitTableau_->describe(out, verbLevel);
875 out <<
" xTilde_ = " << xTilde_ << std::endl;
876 out <<
" stageF_.size() = " << stageF_.size() << std::endl;
877 int numStages = stageF_.size();
878 for (
int i=0; i<numStages; ++i)
879 out <<
" stageF_["<<i<<
"] = " << stageF_[i] << std::endl;
880 out <<
" stageG_.size() = " << stageG_.size() << std::endl;
881 numStages = stageG_.size();
882 for (
int i=0; i<numStages; ++i)
883 out <<
" stageG_["<<i<<
"] = " << stageG_[i] << std::endl;
884 out <<
" stepperRKAppAction_= " << this->stepperRKAppAction_ << std::endl;
885 out <<
" order_ = " << order_ << std::endl;
886 out <<
"--------------------------------" << std::endl;
890template<
class Scalar>
893 out.setOutputToRootOnly(0);
894 bool isValidSetup =
true;
898 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
899 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
900 this->wrapperModel_);
902 if ( wrapperModelPairIMEX->getExplicitModel() == Teuchos::null) {
903 isValidSetup =
false;
904 out <<
"The explicit ModelEvaluator is not set!\n";
907 if ( wrapperModelPairIMEX->getImplicitModel() == Teuchos::null) {
908 isValidSetup =
false;
909 out <<
"The implicit ModelEvaluator is not set!\n";
912 if (this->wrapperModel_ == Teuchos::null) {
913 isValidSetup =
false;
914 out <<
"The wrapper ModelEvaluator is not set!\n";
917 if (this->stepperRKAppAction_ == Teuchos::null) {
918 isValidSetup =
false;
919 out <<
"The AppAction is not set!\n";
922 if ( explicitTableau_ == Teuchos::null ) {
923 isValidSetup =
false;
924 out <<
"The explicit tableau is not set!\n";
927 if ( implicitTableau_ == Teuchos::null ) {
928 isValidSetup =
false;
929 out <<
"The implicit tableau is not set!\n";
936template<
class Scalar>
937Teuchos::RCP<const Teuchos::ParameterList>
940 auto pl = this->getValidParametersBasicImplicit();
941 pl->template set<int>(
"overall order", this->getOrder());
944 explicitStepper->setTableau(
945 explicitTableau_->A(), explicitTableau_->b(), explicitTableau_->c(),
946 explicitTableau_->order(), explicitTableau_->orderMin(),
947 explicitTableau_->orderMax(), explicitTableau_->bstar() );
948 pl->set(
"IMEX-RK Explicit Stepper", *explicitStepper->getValidParameters());
951 implicitStepper->setTableau(
952 implicitTableau_->A(), implicitTableau_->b(), implicitTableau_->c(),
953 implicitTableau_->order(), implicitTableau_->orderMin(),
954 implicitTableau_->orderMax(), implicitTableau_->bstar() );
955 pl->set(
"IMEX-RK Implicit Stepper", *implicitStepper->getValidParameters());
963template<
class Scalar>
964Teuchos::RCP<StepperIMEX_RK<Scalar> >
967 std::string stepperType,
968 Teuchos::RCP<Teuchos::ParameterList> pl)
972 stepper->setStepperImplicitValues(pl);
973 stepper->setTableaus(pl, stepperType);
975 if (model != Teuchos::null) {
976 stepper->setModel(model);
977 stepper->initialize();
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
RK Explicit 3 Stage 3rd order TVD.
General Explicit Runge-Kutta Butcher Tableau.
Time-derivative interface for IMEX RK.
Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void setExplicitTableau(Teuchos::RCP< const RKButcherTableau< Scalar > > explicitTableau)
Set the explicit tableau from tableau.
virtual void initialize()
Initialize during construction and after changing input parameters.
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
virtual void setTableaus(std::string stepperType="", Teuchos::RCP< const RKButcherTableau< Scalar > > explicitTableau=Teuchos::null, Teuchos::RCP< const RKButcherTableau< Scalar > > implicitTableau=Teuchos::null)
Set both the explicit and implicit tableau from ParameterList.
virtual void setImplicitTableau(Teuchos::RCP< const RKButcherTableau< Scalar > > implicitTableau)
Set the implicit tableau from tableau.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Provide a StepperState to the SolutionState. This Stepper does not have any special state data,...
void evalExplicitModel(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &F) const
virtual void setModelPair(const Teuchos::RCP< WrapperModelEvaluatorPairIMEX_Basic< Scalar > > &mePair)
Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set the model.
StepperIMEX_RK(std::string stepperType="IMEX RK SSP2")
Default constructor.
void evalImplicitModelExplicitly(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &G) const
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Application Action for StepperRKBase.
StepperState is a simple class to hold state information about the stepper.
Thyra Base interface for time steppers.
virtual void initialize()
Initialize after construction and changing input parameters.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Teuchos::RCP< StepperIMEX_RK< Scalar > > createStepperIMEX_RK(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, std::string stepperType, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
void validImplicitODE_DAE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate ME supports implicit ODE/DAE evaluation, f(xdot,x,t) [= 0].
void validExplicitODE(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Validate that the model supports explicit ODE evaluation, f(x,t) [=xdot].
@ SOLVE_FOR_X
Solve for x and determine xDot from x.