Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperIMEX_RK_Partition_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_StepperIMEX_RK_Partition_impl_hpp
10#define Tempus_StepperIMEX_RK_Partition_impl_hpp
11
12#include "Thyra_VectorStdOps.hpp"
13
14#include "Tempus_StepperFactory.hpp"
16
17
18namespace Tempus {
19
20
21template<class Scalar>
23{
24 TEUCHOS_TEST_FOR_EXCEPTION(
25 stepperType != "Partitioned IMEX RK 1st order" &&
26 stepperType != "Partitioned IMEX RK SSP2" &&
27 stepperType != "Partitioned IMEX RK ARS 233" &&
28 stepperType != "General Partitioned IMEX RK", std::logic_error,
29 " 'Stepper Type' (='" + stepperType +"')\n"
30 " does not match one of the types for this Stepper:\n"
31 " 'Partitioned IMEX RK 1st order'\n"
32 " 'Partitioned IMEX RK SSP2'\n"
33 " 'Partitioned IMEX RK ARS 233'\n"
34 " 'General Partitioned IMEX RK'\n");
35
36 this->setStepperName( "Partitioned IMEX RK SSP2");
37 this->setStepperType( "Partitioned IMEX RK SSP2");
38 this->setUseFSAL( false);
39 this->setICConsistency( "None");
40 this->setICConsistencyCheck( false);
41 this->setZeroInitialGuess( false);
42
43 this->setStageNumber(-1);
44
45 this->setTableaus(stepperType);
46 this->setAppAction(Teuchos::null);
47 this->setDefaultSolver();
48}
49
50
51template<class Scalar>
53 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
54 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
55 bool useFSAL,
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,
63 Scalar order)
64{
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");
76
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);
84
85 this->setStageNumber(-1);
86
87 if ( stepperType == "General Partitioned IMEX RK" ) {
88 this->setExplicitTableau(explicitTableau);
89 this->setImplicitTableau(implicitTableau);
90 } else {
91 this->setTableaus(stepperType);
92 }
93 this->setAppAction(stepperRKAppAction);
94 this->setSolver(solver);
95
96 if (appModel != Teuchos::null) {
97 this->setModel(appModel);
98 this->initialize();
99 }
100}
101
102
103template<class Scalar>
105 std::string stepperType,
106 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
107 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
108{
109 if (stepperType == "") stepperType = "Partitioned IMEX RK SSP2";
110
111 if (stepperType == "Partitioned IMEX RK 1st order") {
112 {
113 // Explicit Tableau
114 typedef Teuchos::ScalarTraits<Scalar> ST;
115 int NumStages = 2;
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();
121
122 // Fill A:
123 A(0,0) = zero; A(0,1) = zero;
124 A(1,0) = one; A(1,1) = zero;
125
126 // Fill b:
127 b(0) = one; b(1) = zero;
128
129 // Fill c:
130 c(0) = zero; c(1) = one;
131
132 int order = 1;
133
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);
139
140 this->setExplicitTableau(expTableau);
141 }
142 {
143 // Implicit Tableau
144 typedef Teuchos::ScalarTraits<Scalar> ST;
145 int NumStages = 2;
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();
152
153 // Fill A:
154 A(0,0) = zero; A(0,1) = zero;
155 A(1,0) = zero; A(1,1) = one;
156
157 // Fill b:
158 b(0) = zero; b(1) = one;
159
160 // Fill c:
161 c(0) = zero; c(1) = one;
162
163 int order = 1;
164
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);
170
171 this->setImplicitTableau(impTableau);
172 }
173 this->setStepperName("Partitioned IMEX RK 1st order");
174 this->setStepperType("Partitioned IMEX RK 1st order");
175 this->setOrder(1);
176
177 } else if (stepperType == "Partitioned IMEX RK SSP2") {
178 // Explicit Tableau
179 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
180 this->setExplicitTableau(stepperERK->getTableau());
181
182 // Implicit Tableau
183 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
184 stepperSDIRK->setGammaType("2nd Order L-stable");
185 this->setImplicitTableau(stepperSDIRK->getTableau());
186
187 this->setStepperName("Partitioned IMEX RK SSP2");
188 this->setStepperType("Partitioned IMEX RK SSP2");
189 this->setOrder(2);
190
191 } else if (stepperType == "Partitioned IMEX RK ARS 233") {
192 typedef Teuchos::ScalarTraits<Scalar> ST;
193 int NumStages = 3;
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);
201 {
202 // Explicit Tableau
203 // Fill A:
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;
207
208 // Fill b:
209 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
210
211 // Fill c:
212 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
213
214 int order = 2;
215
216 auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
217 "Explicit Tableau - Partitioned IMEX RK ARS 233",
218 A,b,c,order,order,order));
219
220 this->setExplicitTableau(expTableau);
221 }
222 {
223 // Implicit Tableau
224 // Fill A:
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;
228
229 // Fill b:
230 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
231
232 // Fill c:
233 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
234
235 int order = 3;
236
237 auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
238 "Implicit Tableau - Partitioned IMEX RK ARS 233",
239 A,b,c,order,order,order));
240
241 this->setImplicitTableau(impTableau);
242 }
243 this->setStepperName("Partitioned IMEX RK ARS 233");
244 this->setStepperType("Partitioned IMEX RK ARS 233");
245 this->setOrder(3);
246
247 } else if (stepperType == "General Partitioned IMEX RK") {
248
249 if ( explicitTableau == Teuchos::null ) {
250 // Default Explicit Tableau (i.e., Partitioned IMEX RK SSP2)
251 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
252 this->setExplicitTableau(stepperERK->getTableau());
253 } else {
254 this->setExplicitTableau(explicitTableau);
255 }
256
257 if ( implicitTableau == Teuchos::null ) {
258 // Default Implicit Tablea (i.e., Partitioned IMEX RK SSP2)
259 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
260 stepperSDIRK->setGammaType("2nd Order L-stable");
261 this->setImplicitTableau(stepperSDIRK->getTableau());
262 } else {
263 this->setImplicitTableau(implicitTableau);
264 }
265
266 this->setStepperName("General Partitioned IMEX RK");
267 this->setStepperType("General Partitioned IMEX RK");
268 this->setOrder(1);
269
270 } else {
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");
279 }
280
281 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_==Teuchos::null,
282 std::runtime_error,
283 "Error - StepperIMEX_RK_Partition - Explicit tableau is null!");
284 TEUCHOS_TEST_FOR_EXCEPTION(implicitTableau_==Teuchos::null,
285 std::runtime_error,
286 "Error - StepperIMEX_RK_Partition - Implicit tableau is null!");
287 TEUCHOS_TEST_FOR_EXCEPTION(
288 explicitTableau_->numStages()!=implicitTableau_->numStages(),
289 std::runtime_error,
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");
295
296 this->isInitialized_ = false;
297}
298
299
300template<class Scalar>
301void
304 Teuchos::RCP<Teuchos::ParameterList> pl,
305 std::string stepperType)
306{
307 using Teuchos::RCP;
308 if (stepperType == "") {
309 if (pl == Teuchos::null)
310 stepperType = "Partitioned IMEX RK SSP2";
311 else
312 stepperType = pl->get<std::string>("Stepper Type", "Partitioned IMEX RK SSP2");
313 }
314
315 if (stepperType != "General Partitioned IMEX RK") {
316 this->setTableaus(stepperType);
317 } else {
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> > (
327 stepperTemp,true);
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();
331 }
332
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> > (
339 stepperTemp,true);
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();
343 }
344
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");
349
350 this->setTableaus(stepperType, explicitTableau, implicitTableau);
351
352 this->setOrder(pl->get<int>("overall order", 1));
353 }
354 }
355}
356
357
358template<class Scalar>
360 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau)
361{
362 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau->isImplicit() == true,
363 std::logic_error,
364 "Error - Received an implicit Tableau for setExplicitTableau()!\n" <<
365 " Tableau = " << explicitTableau->description() << "\n");
366 explicitTableau_ = explicitTableau;
367
368 this->isInitialized_ = false;
369}
370
371
372template<class Scalar>
374 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
375{
376 TEUCHOS_TEST_FOR_EXCEPTION( implicitTableau->isDIRK() != true,
377 std::logic_error,
378 "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n" <<
379 " Tableau = " << implicitTableau->description() << "\n");
380 implicitTableau_ = implicitTableau;
381
382 this->isInitialized_ = false;
383}
384
385template<class Scalar>
387 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
388{
389 using Teuchos::RCP;
390 using Teuchos::rcp_const_cast;
391 using Teuchos::rcp_dynamic_cast;
392 RCP<Thyra::ModelEvaluator<Scalar> > ncModel =
393 rcp_const_cast<Thyra::ModelEvaluator<Scalar> > (appModel);
394 RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> > modelPairIMEX =
395 rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >(ncModel);
396 TEUCHOS_TEST_FOR_EXCEPTION( modelPairIMEX == Teuchos::null, std::logic_error,
397 "Error - StepperIMEX_RK::setModel() was given a ModelEvaluator that\n"
398 " could not be cast to a WrapperModelEvaluatorPairPartIMEX_Basic!\n"
399 " From: " << appModel << "\n"
400 " To : " << modelPairIMEX << "\n"
401 " Likely have given the wrong ModelEvaluator to this Stepper.\n");
402
403 setModelPair(modelPairIMEX);
404
405 this->isInitialized_ = false;
406}
407
413template<class Scalar>
415 const Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> > &
416 mePairIMEX)
417{
418 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
419 wrapperModelPairIMEX =
420 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
421 (this->wrapperModel_);
422 validExplicitODE (mePairIMEX->getExplicitModel());
423 validImplicitODE_DAE(mePairIMEX->getImplicitModel());
424 wrapperModelPairIMEX = mePairIMEX;
425 wrapperModelPairIMEX->initialize();
426
427 this->wrapperModel_ = wrapperModelPairIMEX;
428
429 this->isInitialized_ = false;
430}
431
437template<class Scalar>
439 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
440 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel)
441{
442 validExplicitODE (explicitModel);
443 validImplicitODE_DAE(implicitModel);
444 this->wrapperModel_ = Teuchos::rcp(
445 new WrapperModelEvaluatorPairPartIMEX_Basic<Scalar>(
446 explicitModel, implicitModel));
447
448 this->isInitialized_ = false;
449}
450
451
452template<class Scalar>
454{
455 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
456 wrapperModelPairIMEX =
457 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
458 (this->wrapperModel_);
459 TEUCHOS_TEST_FOR_EXCEPTION( wrapperModelPairIMEX == Teuchos::null,
460 std::logic_error,
461 "Error - Can not cast the wrapper Model Evaluator to a IMEX Model Pair."
462 "StepperIMEX_RK_Partition::initialize()\n");
463
464 // Initialize the stage vectors
465 const int numStages = explicitTableau_->numStages();
466 stageF_.resize(numStages);
467 stageGx_.resize(numStages);
468 for(int i=0; i < numStages; i++) {
469 stageF_[i] = Thyra::createMember(wrapperModelPairIMEX->
470 getExplicitModel()->get_f_space());
471 stageGx_[i] = Thyra::createMember(wrapperModelPairIMEX->
472 getImplicitModel()->get_f_space());
473 assign(stageF_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
474 assign(stageGx_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
475 }
476
477 xTilde_ = Thyra::createMember(wrapperModelPairIMEX->
478 getImplicitModel()->get_x_space());
479 assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
480
482}
483
484
485template<class Scalar>
487 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
488{
489 using Teuchos::RCP;
490
491 int numStates = solutionHistory->getNumStates();
492
493 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
494 "Error - setInitialConditions() needs at least one SolutionState\n"
495 " to set the initial condition. Number of States = " << numStates);
496
497 if (numStates > 1) {
498 RCP<Teuchos::FancyOStream> out = this->getOStream();
499 Teuchos::OSTab ostab(out,1,"StepperIMEX_RK::setInitialConditions()");
500 *out << "Warning -- SolutionHistory has more than one state!\n"
501 << "Setting the initial conditions on the currentState.\n"<<std::endl;
502 }
503
504 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
505 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
506
507 // Use x from inArgs as ICs, if needed.
508 auto inArgs = this->wrapperModel_->getNominalValues();
509 if (x == Teuchos::null) {
510 TEUCHOS_TEST_FOR_EXCEPTION( (x == Teuchos::null) &&
511 (inArgs.get_x() == Teuchos::null), std::logic_error,
512 "Error - setInitialConditions() needs the ICs from the SolutionHistory\n"
513 " or getNominalValues()!\n");
514
515 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
516 initialState->setX(x);
517 }
518
519
520 // Perform IC Consistency
521 std::string icConsistency = this->getICConsistency();
522 TEUCHOS_TEST_FOR_EXCEPTION(icConsistency != "None", std::logic_error,
523 "Error - setInitialConditions() requested a consistency of '"
524 << icConsistency << "'.\n"
525 " But only 'None' is available for IMEX-RK!\n");
526
527 TEUCHOS_TEST_FOR_EXCEPTION( this->getUseFSAL(), std::logic_error,
528 "Error - The First-Same-As-Last (FSAL) principle is not "
529 << "available for IMEX-RK. Set useFSAL=false.\n");
530}
531
532
533template <typename Scalar>
535 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & X,
536 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & Y,
537 Scalar time, Scalar stepSize, Scalar stageNumber,
538 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & G) const
539{
540 typedef Thyra::ModelEvaluatorBase MEB;
541 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
542 wrapperModelPairIMEX =
543 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
544 (this->wrapperModel_);
545 MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->getInArgs();
546 inArgs.set_x(X);
547 inArgs.set_p(wrapperModelPairIMEX->getParameterIndex(), Y);
548 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
549 if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
550 if (inArgs.supports(MEB::IN_ARG_stage_number))
551 inArgs.set_stage_number(stageNumber);
552
553 // For model evaluators whose state function f(x, x_dot, t) describes
554 // an implicit ODE, and which accept an optional x_dot input argument,
555 // make sure the latter is set to null in order to request the evaluation
556 // of a state function corresponding to the explicit ODE formulation
557 // x_dot = f(x, t)
558 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
559
560 MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->getOutArgs();
561 outArgs.set_f(G);
562
563 wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs,outArgs);
564 Thyra::Vt_S(G.ptr(), -1.0);
565}
566
567
568template <typename Scalar>
570 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & Z,
571 Scalar time, Scalar stepSize, Scalar stageNumber,
572 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & F) const
573{
574 typedef Thyra::ModelEvaluatorBase MEB;
575
576 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
577 wrapperModelPairIMEX =
578 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
579 (this->wrapperModel_);
580 MEB::InArgs<Scalar> inArgs =
581 wrapperModelPairIMEX->getExplicitModel()->createInArgs();
582 inArgs.set_x(Z);
583 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
584 if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
585 if (inArgs.supports(MEB::IN_ARG_stage_number))
586 inArgs.set_stage_number(stageNumber);
587
588 // For model evaluators whose state function f(x, x_dot, t) describes
589 // an implicit ODE, and which accept an optional x_dot input argument,
590 // make sure the latter is set to null in order to request the evaluation
591 // of a state function corresponding to the explicit ODE formulation
592 // x_dot = f(x, t)
593 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
594
595 MEB::OutArgs<Scalar> outArgs =
596 wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
597 outArgs.set_f(F);
598
599 wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
600 Thyra::Vt_S(F.ptr(), -1.0);
601}
602
603
604template<class Scalar>
606 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
607{
608 this->checkInitialized();
609
610 using Teuchos::RCP;
611 using Teuchos::SerialDenseMatrix;
612 using Teuchos::SerialDenseVector;
613
614 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperIMEX_RK_Partition::takeStep()");
615 {
616 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
617 std::logic_error,
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");
623
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();
628
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();
636
637 Teuchos::RCP<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
638 wrapperModelPairIMEX =
639 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairPartIMEX_Basic<Scalar> >
640 (this->wrapperModel_);
641
642 bool pass = true;
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());
648
649 RCP<StepperIMEX_RK_Partition<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
650 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
652
653 // Compute stage solutions
654 for (int i = 0; i < numStages; ++i) {
655 this->setStageNumber(i);
656
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);
669 }
670 if (A (i,j) != Teuchos::ScalarTraits<Scalar>::zero())
671 Thyra::Vp_StV(xTilde_.ptr(), -dt*A (i,j), *(stageGx_[j]));
672 }
673
674 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
676
677 Scalar ts = time + c(i)*dt;
678 Scalar tHats = time + cHat(i)*dt;
679 if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
680 // Explicit stage for the ImplicitODE_DAE
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) {
685 // stageGx_[i] is not needed.
686 assign(stageGx_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
687 } else {
688 Thyra::assign(stageX.ptr(), *xTilde_);
689 evalImplicitModelExplicitly(stageX, stageY, ts, dt, i, stageGx_[i]);
690 }
691 } else {
692 // Implicit stage for the ImplicitODE_DAE
693 const Scalar alpha = Scalar(1.0)/(dt*A(i,i));
694 const Scalar beta = Scalar(1.0);
695
696 // Setup TimeDerivative
697 RCP<TimeDerivative<Scalar> > timeDer =
698 Teuchos::rcp(new StepperIMEX_RKPartTimeDerivative<Scalar>(
699 alpha, xTilde_.getConst()));
700
701 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
702 timeDer, dt, alpha, beta, SOLVE_FOR_X, i));
703
704 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
706
707 wrapperModelPairIMEX->setUseImplicitModel(true);
708 this->solver_->setModel(wrapperModelPairIMEX);
709
710 Thyra::SolveStatus<Scalar> sStatus;
711 if (wrapperModelPairIMEX->getParameterIndex() >= 0)
712 sStatus = this->solveImplicitODE(stageX, stageGx_[i], ts, p,
713 stageY, wrapperModelPairIMEX->getParameterIndex());
714 else
715 sStatus = this->solveImplicitODE(stageX, stageGx_[i], ts, p);
716
717 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass = false;
718
719 wrapperModelPairIMEX->setUseImplicitModel(false);
720
721 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
723
724 // Update contributions to stage values
725 Thyra::V_StVpStV(stageGx_[i].ptr(), -alpha, *stageX, alpha, *xTilde_);
726 }
727
728 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
730 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
731 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
733 }
734
735 // Sum for solution: y_n = y_n-1 - dt*Sum{ bHat(i)*fy(i) }
736 // Sum for solution: x_n = x_n-1 - dt*Sum{ bHat(i)*fx(i) + b(i)*gx(i) }
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]));
745 }
746
747 if (pass == true) workingState->setSolutionStatus(Status::PASSED);
748 else workingState->setSolutionStatus(Status::FAILED);
749 workingState->setOrder(this->getOrder());
750 workingState->computeNorms(currentState);
751 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
753 }
754 // reset the stage number
755 this->setStageNumber(-1);
756 return;
757}
758
765template<class Scalar>
766Teuchos::RCP<Tempus::StepperState<Scalar> >
769{
770 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
771 rcp(new StepperState<Scalar>(this->getStepperType()));
772 return stepperState;
773}
774
775
776template<class Scalar>
778 Teuchos::FancyOStream &out,
779 const Teuchos::EVerbosityLevel verbLevel) const
780{
781 out.setOutputToRootOnly(0);
782 out << std::endl;
783 Stepper<Scalar>::describe(out, verbLevel);
784 StepperImplicit<Scalar>::describe(out, verbLevel);
785
786 out << "--- StepperIMEX_RK_Partition ---\n";
787 out << " explicitTableau_ = " << explicitTableau_ << std::endl;
788 if (verbLevel == Teuchos::VERB_HIGH)
789 explicitTableau_->describe(out, verbLevel);
790 out << " implicitTableau_ = " << implicitTableau_ << std::endl;
791 if (verbLevel == Teuchos::VERB_HIGH)
792 implicitTableau_->describe(out, verbLevel);
793 out << " xTilde_ = " << xTilde_ << std::endl;
794 out << " stageF_.size() = " << stageF_.size() << std::endl;
795 int numStages = stageF_.size();
796 for (int i=0; i<numStages; ++i)
797 out << " stageF_["<<i<<"] = " << stageF_[i] << std::endl;
798 out << " stageGx_.size() = " << stageGx_.size() << std::endl;
799 numStages = stageGx_.size();
800 for (int i=0; i<numStages; ++i)
801 out << " stageGx_["<<i<<"] = " << stageGx_[i] << std::endl;
802 out << " stepperRKAppAction_= " << this->stepperRKAppAction_ << std::endl;
803 out << " order_ = " << order_ << std::endl;
804 out << "--------------------------------" << std::endl;
805}
806
807
808template<class Scalar>
809bool StepperIMEX_RK_Partition<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
810{
811 out.setOutputToRootOnly(0);
812
813 bool isValidSetup = true;
814
815 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
816
817 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
818 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
819 this->wrapperModel_);
820
821 if ( wrapperModelPairIMEX->getExplicitModel() == Teuchos::null) {
822 isValidSetup = false;
823 out << "The explicit ModelEvaluator is not set!\n";
824 }
825
826 if ( wrapperModelPairIMEX->getImplicitModel() == Teuchos::null) {
827 isValidSetup = false;
828 out << "The implicit ModelEvaluator is not set!\n";
829 }
830
831 if (this->wrapperModel_ == Teuchos::null) {
832 isValidSetup = false;
833 out << "The wrapper ModelEvaluator is not set!\n";
834 }
835
836 if (this->solver_ == Teuchos::null) {
837 isValidSetup = false;
838 out << "The solver is not set!\n";
839 }
840
841 if (this->stepperRKAppAction_ == Teuchos::null) {
842 isValidSetup = false;
843 out << "The AppAction is not set!\n";
844 }
845
846 if ( explicitTableau_ == Teuchos::null ) {
847 isValidSetup = false;
848 out << "The explicit tableau is not set!\n";
849 }
850
851 if ( implicitTableau_ == Teuchos::null ) {
852 isValidSetup = false;
853 out << "The implicit tableau is not set!\n";
854 }
855
856 return isValidSetup;
857}
858
859
860template<class Scalar>
861Teuchos::RCP<const Teuchos::ParameterList>
863{
864 auto pl = this->getValidParametersBasicImplicit();
865 pl->template set<int>("overall order", this->getOrder());
866
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());
873
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());
880
881 return pl;
882}
883
884
885// Nonmember constructor - ModelEvaluator and ParameterList
886// ------------------------------------------------------------------------
887template<class Scalar>
888Teuchos::RCP<StepperIMEX_RK_Partition<Scalar> >
890 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
891 std::string stepperType,
892 Teuchos::RCP<Teuchos::ParameterList> pl)
893{
894 auto stepper = Teuchos::rcp(new StepperIMEX_RK_Partition<Scalar>(stepperType));
895 stepper->setStepperImplicitValues(pl);
896 stepper->setTableausPartition(pl, stepperType);
897
898 if (model != Teuchos::null) {
899 stepper->setModel(model);
900 stepper->initialize();
901 }
902
903 return stepper;
904}
905
906
907} // namespace Tempus
908#endif // Tempus_StepperIMEX_RK_Partition_impl_hpp
virtual void setModelPair(const Teuchos::RCP< WrapperModelEvaluatorPairPartIMEX_Basic< Scalar > > &modelPair)
Create WrapperModelPairIMEX from user-supplied ModelEvaluator pair.
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.
StepperIMEX_RK_Partition(std::string stepperType="Partitioned IMEX RK SSP2")
Default constructor.
virtual Teuchos::RCP< Tempus::StepperState< Scalar > > getDefaultStepperState()
Provide a StepperState to the SolutionState. This Stepper does not have any special state data,...
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void setExplicitTableau(Teuchos::RCP< const RKButcherTableau< Scalar > > explicitTableau)
Set the explicit tableau from tableau.
virtual void setInitialConditions(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Set the initial conditions and make them consistent.
virtual void setModel(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &appModel)
Set the model.
virtual void setTableausPartition(Teuchos::RCP< Teuchos::ParameterList > pl, std::string stepperType)
virtual bool isValidSetup(Teuchos::FancyOStream &out) const
void evalImplicitModelExplicitly(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &X, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &Y, Scalar time, Scalar stepSize, Scalar stageNumber, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &G) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual void initialize()
Initialize during construction and after changing input parameters.
virtual void setImplicitTableau(Teuchos::RCP< const RKButcherTableau< Scalar > > implicitTableau)
Set the implicit tableau from tableau.
virtual void takeStep(const Teuchos::RCP< SolutionHistory< Scalar > > &solutionHistory)
Take the specified timestep, dt, and return true if successful.
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 describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
Application Action for StepperRKBase.
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
Teuchos::RCP< StepperIMEX_RK_Partition< Scalar > > createStepperIMEX_RK_Partition(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.