Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_StepperIMEX_RK_impl.hpp
Go to the documentation of this file.
1//@HEADER
2// *****************************************************************************
3// Tempus: Time Integration and Sensitivity Analysis Package
4//
5// Copyright 2017 NTESS and the Tempus contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8//@HEADER
9
10#ifndef Tempus_StepperIMEX_RK_impl_hpp
11#define Tempus_StepperIMEX_RK_impl_hpp
12
13#include "Thyra_VectorStdOps.hpp"
14
15#include "Tempus_StepperFactory.hpp"
17
18namespace Tempus {
19
20template <class Scalar>
22{
23 TEUCHOS_TEST_FOR_EXCEPTION(
24 stepperType != "IMEX RK 1st order" && stepperType != "SSP1_111" &&
25 stepperType != "IMEX RK SSP2" && stepperType != "IMEX RK SSP3" &&
26 stepperType != "SSP3_332" && stepperType != "SSP2_222" &&
27 stepperType != "SSP2_222_L" && stepperType != "SSP2_222_A" &&
28 stepperType != "IMEX RK ARS 233" && stepperType != "ARS 233" &&
29 stepperType != "General IMEX RK",
30 std::logic_error,
31 " 'Stepper Type' (='"
32 << stepperType
33 << "')\n"
34 << " does not match one of the types for this Stepper:\n"
35 << " 'IMEX RK 1st order'\n"
36 << " 'SSP1_111'\n"
37 << " 'IMEX RK SSP2'\n"
38 << " 'IMEX RK SSP3'\n"
39 << " 'SSP3_332'\n"
40 << " 'SSP2_222'\n"
41 << " 'SSP2_222_L'\n"
42 << " 'SSP2_222_A'\n"
43 << " 'IMEX RK ARS 233'\n"
44 << " 'ARS 233'\n"
45 << " 'General IMEX RK'\n");
46
47 this->setStepperName(stepperType);
48 this->setStepperType(stepperType);
49 this->setUseFSAL(false);
50 this->setICConsistency("None");
51 this->setICConsistencyCheck(false);
52 this->setZeroInitialGuess(false);
53
54 this->setStageNumber(-1);
55
56 this->setTableaus(stepperType);
57 this->setAppAction(Teuchos::null);
58 this->setDefaultSolver();
59}
60
61template <class Scalar>
63 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
64 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
65 bool useFSAL, std::string ICConsistency, bool ICConsistencyCheck,
66 bool zeroInitialGuess,
67 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
68 std::string stepperType,
69 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
70 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau, Scalar order)
71{
72 TEUCHOS_TEST_FOR_EXCEPTION(
73 stepperType != "IMEX RK 1st order" && stepperType != "SSP1_111" &&
74 stepperType != "IMEX RK SSP2" && stepperType != "IMEX RK SSP3" &&
75 stepperType != "SSP3_332" && stepperType != "SSP2_222" &&
76 stepperType != "SSP2_222_L" && stepperType != "SSP2_222_A" &&
77 stepperType != "IMEX RK ARS 233" && stepperType != "ARS 233" &&
78 stepperType != "General IMEX RK",
79 std::logic_error,
80 " 'Stepper Type' (='"
81 << stepperType
82 << "')\n"
83 << " does not match one of the types for this Stepper:\n"
84 << " 'IMEX RK 1st order'\n"
85 << " 'SSP1_111'\n"
86 << " 'IMEX RK SSP2'\n"
87 << " 'IMEX RK SSP3'\n"
88 << " 'SSP3_332'\n"
89 << " 'SSP2_222'\n"
90 << " 'SSP2_222_L'\n"
91 << " 'SSP2_222_A'\n"
92 << " 'IMEX RK ARS 233'\n"
93 << " 'ARS 233'\n"
94 << " 'General IMEX RK'\n");
95
96 this->setStepperName(stepperType);
97 this->setStepperType(stepperType);
98 this->setUseFSAL(useFSAL);
99 this->setICConsistency(ICConsistency);
100 this->setICConsistencyCheck(ICConsistencyCheck);
101 this->setZeroInitialGuess(zeroInitialGuess);
102
103 this->setStageNumber(-1);
104
105 if (stepperType == "General IMEX RK") {
106 this->setExplicitTableau(explicitTableau);
107 this->setImplicitTableau(implicitTableau);
108 }
109 else {
110 this->setTableaus(stepperType);
111 }
112 this->setOrder(order);
113 this->setAppAction(stepperRKAppAction);
114 this->setSolver(solver);
115
116 if (appModel != Teuchos::null) {
117 this->setModel(appModel);
118 this->initialize();
119 }
120}
121
122template <class Scalar>
124 std::string stepperType,
125 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
126 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
127{
128 if (stepperType == "") stepperType = "IMEX RK SSP2";
129
130 if (stepperType == "IMEX RK 1st order") {
131 {
132 // Explicit Tableau
133 typedef Teuchos::ScalarTraits<Scalar> ST;
134 int NumStages = 2;
135 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
136 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
137 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
138 const Scalar one = ST::one();
139 const Scalar zero = ST::zero();
140
141 // Fill A:
142 A(0, 0) = zero;
143 A(0, 1) = zero;
144 A(1, 0) = one;
145 A(1, 1) = zero;
146
147 // Fill b:
148 b(0) = one;
149 b(1) = zero;
150
151 // Fill c:
152 c(0) = zero;
153 c(1) = one;
154
155 int order = 1;
156
157 auto expTableau = Teuchos::rcp(
158 new RKButcherTableau<Scalar>("Explicit Tableau - IMEX RK 1st order",
159 A, b, c, order, order, order));
160 expTableau->setTVD(true);
161 expTableau->setTVDCoeff(2.0);
162
163 this->setExplicitTableau(expTableau);
164 }
165 {
166 // Implicit Tableau
167 typedef Teuchos::ScalarTraits<Scalar> ST;
168 int NumStages = 2;
169 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
170 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
171 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
172 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
173 const Scalar one = ST::one();
174 const Scalar zero = ST::zero();
175
176 // Fill A:
177 A(0, 0) = zero;
178 A(0, 1) = zero;
179 A(1, 0) = zero;
180 A(1, 1) = one;
181
182 // Fill b:
183 b(0) = zero;
184 b(1) = one;
185
186 // Fill c:
187 c(0) = zero;
188 c(1) = one;
189
190 int order = 1;
191
192 auto impTableau = Teuchos::rcp(
193 new RKButcherTableau<Scalar>("Implicit Tableau - IMEX RK 1st order",
194 A, b, c, order, order, order));
195 impTableau->setTVD(true);
196 impTableau->setTVDCoeff(sspcoef);
197
198 this->setImplicitTableau(impTableau);
199 }
200 this->setStepperName("IMEX RK 1st order");
201 this->setStepperType("IMEX RK 1st order");
202 this->setOrder(1);
203 }
204 else if (stepperType == "SSP1_111") {
205 {
206 // Explicit Tableau
207 typedef Teuchos::ScalarTraits<Scalar> ST;
208 const int NumStages = 1;
209 const int order = 1;
210 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
211 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
212 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
213 const Scalar one = ST::one();
214 const Scalar zero = ST::zero();
215
216 // Fill A:
217 A(0, 0) = zero;
218
219 // Fill b:
220 b(0) = one;
221
222 // Fill c:
223 c(0) = zero;
224
225 auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
226 "Explicit Tableau - SSP1_111", A, b, c, order, order, order));
227 expTableau->setTVD(true);
228 expTableau->setTVDCoeff(1.0);
229
230 this->setExplicitTableau(expTableau);
231 }
232 {
233 // Implicit Tableau
234 typedef Teuchos::ScalarTraits<Scalar> ST;
235 const int NumStages = 1;
236 const int order = 1;
237 const Scalar sspcoef = std::numeric_limits<Scalar>::max();
238 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
239 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
240 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
241 const Scalar one = ST::one();
242
243 // Fill A:
244 A(0, 0) = one;
245
246 // Fill b:
247 b(0) = one;
248
249 // Fill c:
250 c(0) = one;
251
252 auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
253 "Implicit Tableau - SSP1_111", A, b, c, order, order, order));
254 impTableau->setTVD(true);
255 impTableau->setTVDCoeff(sspcoef);
256
257 this->setImplicitTableau(impTableau);
258 }
259 this->setStepperName("SSP1_111");
260 this->setStepperType("SSP1_111");
261 this->setOrder(1);
262 }
263 else if (stepperType == "IMEX RK SSP2" || stepperType == "SSP2_222_L") {
264 // Explicit Tableau
265 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
266 this->setExplicitTableau(stepperERK->getTableau());
267
268 // Implicit Tableau
269 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
270 stepperSDIRK->setGammaType("2nd Order L-stable");
271 this->setImplicitTableau(stepperSDIRK->getTableau());
272
273 this->setStepperName("IMEX RK SSP2");
274 this->setStepperType("IMEX RK SSP2");
275 this->setOrder(2);
276 }
277 else if (stepperType == "SSP2_222" || stepperType == "SSP2_222_A") {
278 // Explicit Tableau
279 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
280 this->setExplicitTableau(stepperERK->getTableau());
281
282 // Implicit Tableau
283 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
284 stepperSDIRK->setGammaType("gamma");
285 stepperSDIRK->setGamma(0.5);
286 this->setImplicitTableau(stepperSDIRK->getTableau());
287
288 this->setStepperName("SSP2_222");
289 this->setStepperType("SSP2_222");
290 this->setOrder(2);
291 }
292 else if (stepperType == "IMEX RK SSP3" || stepperType == "SSP3_332") {
293 // Explicit Tableau
294 auto stepperERK = Teuchos::rcp(new StepperERK_3Stage3rdOrderTVD<Scalar>());
295 this->setExplicitTableau(stepperERK->getTableau());
296
297 // Implicit Tableau
298 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_3Stage2ndOrder<Scalar>());
299 this->setImplicitTableau(stepperSDIRK->getTableau());
300
301 this->setStepperName("IMEX RK SSP3");
302 this->setStepperType("IMEX RK SSP3");
303 this->setOrder(2);
304 }
305 else if (stepperType == "IMEX RK ARS 233" || stepperType == "ARS 233") {
306 typedef Teuchos::ScalarTraits<Scalar> ST;
307 int NumStages = 3;
308 Teuchos::SerialDenseMatrix<int, Scalar> A(NumStages, NumStages);
309 Teuchos::SerialDenseVector<int, Scalar> b(NumStages);
310 Teuchos::SerialDenseVector<int, Scalar> c(NumStages);
311 const Scalar one = ST::one();
312 const Scalar zero = ST::zero();
313 const Scalar onehalf = ST::one() / (2 * ST::one());
314 const Scalar gamma = (3 * one + ST::squareroot(3 * one)) / (6 * one);
315 {
316 // Explicit Tableau
317 // Fill A:
318 A(0, 0) = zero;
319 A(0, 1) = zero;
320 A(0, 2) = zero;
321 A(1, 0) = gamma;
322 A(1, 1) = zero;
323 A(1, 2) = zero;
324 A(2, 0) = (gamma - 1.0);
325 A(2, 1) = (2.0 - 2.0 * gamma);
326 A(2, 2) = zero;
327
328 // Fill b:
329 b(0) = zero;
330 b(1) = onehalf;
331 b(2) = onehalf;
332
333 // Fill c:
334 c(0) = zero;
335 c(1) = gamma;
336 c(2) = one - gamma;
337
338 int order = 2;
339
340 auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
341 "Partition IMEX-RK Explicit Stepper", A, b, c, order, order, order));
342
343 this->setExplicitTableau(expTableau);
344 }
345 {
346 // Implicit Tableau
347 // Fill A:
348 A(0, 0) = zero;
349 A(0, 1) = zero;
350 A(0, 2) = zero;
351 A(1, 0) = zero;
352 A(1, 1) = gamma;
353 A(1, 2) = zero;
354 A(2, 0) = zero;
355 A(2, 1) = (1.0 - 2.0 * gamma);
356 A(2, 2) = gamma;
357
358 // Fill b:
359 b(0) = zero;
360 b(1) = onehalf;
361 b(2) = onehalf;
362
363 // Fill c:
364 c(0) = zero;
365 c(1) = gamma;
366 c(2) = one - gamma;
367
368 int order = 3;
369
370 auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
371 "Partition IMEX-RK Implicit Stepper", A, b, c, order, order, order));
372
373 this->setImplicitTableau(impTableau);
374 }
375 this->setStepperName("IMEX RK ARS 233");
376 this->setStepperType("IMEX RK ARS 233");
377 this->setOrder(3);
378 }
379 else if (stepperType == "General IMEX RK") {
380 if (explicitTableau == Teuchos::null) {
381 // Default Explicit Tableau (i.e., IMEX RK SSP2)
382 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
383 this->setExplicitTableau(stepperERK->getTableau());
384 }
385 else {
386 this->setExplicitTableau(explicitTableau);
387 }
388
389 if (explicitTableau == Teuchos::null) {
390 // Default Implicit Tableau (i.e., IMEX RK SSP2)
391 auto stepperSDIRK =
392 Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
393 stepperSDIRK->setGammaType("2nd Order L-stable");
394 this->setImplicitTableau(stepperSDIRK->getTableau());
395 }
396 else {
397 this->setImplicitTableau(implicitTableau);
398 }
399
400 this->setStepperName("General IMEX RK");
401 this->setStepperType("General IMEX RK");
402 this->setOrder(1);
403 }
404 else {
405 TEUCHOS_TEST_FOR_EXCEPTION(
406 true, std::logic_error,
407 "Error - Not a valid StepperIMEX_RK type! Stepper Type = "
408 << stepperType << "\n"
409 << " Current valid types are: \n"
410 << " 'IMEX RK 1st order\n"
411 << " 'SSP1_111\n"
412 << " 'IMEX RK SSP2' ('SSP2_222_L')\n"
413 << " 'SSP2_222' ('SSP2_222_A')\n"
414 << " 'IMEX RK SSP3' ('SSP3_332')\n"
415 << " 'IMEX RK ARS 233' ('ARS 233')\n"
416 << " 'General IMEX RK'\n");
417 }
418
419 TEUCHOS_TEST_FOR_EXCEPTION(
420 explicitTableau_ == Teuchos::null, std::runtime_error,
421 "Error - StepperIMEX_RK - Explicit tableau is null!");
422 TEUCHOS_TEST_FOR_EXCEPTION(
423 implicitTableau_ == Teuchos::null, std::runtime_error,
424 "Error - StepperIMEX_RK - Implicit tableau is null!");
425 TEUCHOS_TEST_FOR_EXCEPTION(
426 explicitTableau_->numStages() != implicitTableau_->numStages(),
427 std::runtime_error,
428 "Error - StepperIMEX_RK - Number of stages do not match!\n"
429 << " Explicit tableau = " << explicitTableau_->description() << "\n"
430 << " number of stages = " << explicitTableau_->numStages() << "\n"
431 << " Implicit tableau = " << implicitTableau_->description() << "\n"
432 << " number of stages = " << implicitTableau_->numStages() << "\n");
433
434 this->isInitialized_ = false;
435}
436
437template <class Scalar>
439 Teuchos::RCP<Teuchos::ParameterList> pl, std::string stepperType)
440{
441 using Teuchos::RCP;
442 if (stepperType == "") {
443 if (pl == Teuchos::null)
444 stepperType = "IMEX RK SSP2";
445 else
446 stepperType = pl->get<std::string>("Stepper Type", "IMEX RK SSP2");
447 }
448
449 if (stepperType != "General IMEX RK") {
450 this->setTableaus(stepperType);
451 }
452 else {
453 if (pl != Teuchos::null) {
454 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau;
455 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau;
456 if (pl->isSublist("IMEX-RK Explicit Stepper")) {
457 RCP<Teuchos::ParameterList> explicitPL =
458 Teuchos::rcp(new Teuchos::ParameterList(
459 pl->sublist("IMEX-RK Explicit Stepper")));
460 auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
461 auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
462 auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> >(
463 stepperTemp, true);
464 TEUCHOS_TEST_FOR_EXCEPTION(
465 stepperERK == Teuchos::null, std::logic_error,
466 "Error - The explicit component of a general IMEX RK stepper was "
467 "not specified as an ExplicitRK stepper");
468 explicitTableau = stepperERK->getTableau();
469 }
470
471 if (pl->isSublist("IMEX-RK Implicit Stepper")) {
472 RCP<Teuchos::ParameterList> implicitPL =
473 Teuchos::rcp(new Teuchos::ParameterList(
474 pl->sublist("IMEX-RK Implicit Stepper")));
475 auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
476 auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
477 auto stepperDIRK =
478 Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> >(stepperTemp, true);
479 TEUCHOS_TEST_FOR_EXCEPTION(
480 stepperDIRK == Teuchos::null, std::logic_error,
481 "Error - The implicit component of a general IMEX RK stepper was "
482 "not specified as an DIRK stepper");
483 implicitTableau = stepperDIRK->getTableau();
484 }
485
486 TEUCHOS_TEST_FOR_EXCEPTION(
487 !(explicitTableau != Teuchos::null &&
488 implicitTableau != Teuchos::null),
489 std::logic_error,
490 "Error - A parameter list was used to setup a general IMEX RK "
491 "stepper, but did not "
492 "specify both an explicit and an implicit tableau!\n");
493
494 this->setTableaus(stepperType, explicitTableau, implicitTableau);
495 this->setOrder(pl->get<int>("overall order", 1));
496 }
497 }
498}
499
500template <class Scalar>
502 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau)
503{
504 TEUCHOS_TEST_FOR_EXCEPTION(
505 explicitTableau->isImplicit() == true, std::logic_error,
506 "Error - Received an implicit Tableau for setExplicitTableau()!\n"
507 << " Tableau = " << explicitTableau->description() << "\n");
508 explicitTableau_ = explicitTableau;
509
510 this->isInitialized_ = false;
511}
512
513template <class Scalar>
515 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
516{
517 TEUCHOS_TEST_FOR_EXCEPTION(
518 implicitTableau->isDIRK() != true, std::logic_error,
519 "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n"
520 << " Tableau = " << implicitTableau->description() << "\n");
521 implicitTableau_ = implicitTableau;
522
523 this->isInitialized_ = false;
524}
525
526template <class Scalar>
528 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
529{
530 using Teuchos::RCP;
531 using Teuchos::rcp_const_cast;
532 using Teuchos::rcp_dynamic_cast;
533 RCP<Thyra::ModelEvaluator<Scalar> > ncModel =
534 rcp_const_cast<Thyra::ModelEvaluator<Scalar> >(appModel);
535 RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > modelPairIMEX =
536 rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX_Basic<Scalar> >(ncModel);
537 TEUCHOS_TEST_FOR_EXCEPTION(
538 modelPairIMEX == Teuchos::null, std::logic_error,
539 "Error - StepperIMEX_RK::setModel() was given a ModelEvaluator that\n"
540 << " could not be cast to a WrapperModelEvaluatorPairIMEX_Basic!\n"
541 << " From: " << appModel << "\n"
542 << " To : " << modelPairIMEX << "\n"
543 << " Likely have given the wrong ModelEvaluator to this Stepper.\n");
544
545 setModelPair(modelPairIMEX);
546
547 TEUCHOS_TEST_FOR_EXCEPTION(this->solver_ == Teuchos::null, std::logic_error,
548 "Error - Solver is not set!\n");
549 this->solver_->setModel(modelPairIMEX);
550
551 this->isInitialized_ = false;
552}
553
559template <class Scalar>
562 modelPairIMEX)
563{
564 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
565 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
566 this->wrapperModel_);
567 validExplicitODE(modelPairIMEX->getExplicitModel());
568 validImplicitODE_DAE(modelPairIMEX->getImplicitModel());
569 wrapperModelPairIMEX = modelPairIMEX;
570 wrapperModelPairIMEX->initialize();
571 int expXDim = wrapperModelPairIMEX->getExplicitModel()->get_x_space()->dim();
572 int impXDim = wrapperModelPairIMEX->getImplicitModel()->get_x_space()->dim();
573 TEUCHOS_TEST_FOR_EXCEPTION(
574 expXDim != impXDim, std::logic_error,
575 "Error - \n"
576 << " Explicit and Implicit x vectors are incompatible!\n"
577 << " Explicit vector dim = " << expXDim << "\n"
578 << " Implicit vector dim = " << impXDim << "\n");
579
580 this->wrapperModel_ = wrapperModelPairIMEX;
581
582 this->isInitialized_ = false;
583}
584
590template <class Scalar>
592 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
593 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel)
594{
595 validExplicitODE(explicitModel);
596 validImplicitODE_DAE(implicitModel);
597 this->wrapperModel_ =
599 explicitModel, implicitModel));
600
601 this->isInitialized_ = false;
602}
603
604template <class Scalar>
606{
607 TEUCHOS_TEST_FOR_EXCEPTION(
608 this->wrapperModel_ == Teuchos::null, std::logic_error,
609 "Error - Need to set the model, setModel(), before calling "
610 "StepperIMEX_RK::initialize()\n");
611
612 // Initialize the stage vectors
613 const int numStages = explicitTableau_->numStages();
614 stageF_.resize(numStages);
615 stageG_.resize(numStages);
616 for (int i = 0; i < numStages; i++) {
617 stageF_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
618 stageG_[i] = Thyra::createMember(this->wrapperModel_->get_f_space());
619 assign(stageF_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
620 assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
621 }
622
623 xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
624 assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
625
627}
628
629template <class Scalar>
631 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
632{
633 using Teuchos::RCP;
634
635 int numStates = solutionHistory->getNumStates();
636
637 TEUCHOS_TEST_FOR_EXCEPTION(
638 numStates < 1, std::logic_error,
639 "Error - setInitialConditions() needs at least one SolutionState\n"
640 " to set the initial condition. Number of States = "
641 << numStates);
642
643 if (numStates > 1) {
644 RCP<Teuchos::FancyOStream> out = this->getOStream();
645 Teuchos::OSTab ostab(out, 1, "StepperIMEX_RK::setInitialConditions()");
646 *out << "Warning -- SolutionHistory has more than one state!\n"
647 << "Setting the initial conditions on the currentState.\n"
648 << std::endl;
649 }
650
651 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
652 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
653
654 // Use x from inArgs as ICs, if needed.
655 auto inArgs = this->wrapperModel_->getNominalValues();
656 if (x == Teuchos::null) {
657 TEUCHOS_TEST_FOR_EXCEPTION(
658 (x == Teuchos::null) && (inArgs.get_x() == Teuchos::null),
659 std::logic_error,
660 "Error - setInitialConditions() needs the ICs from the "
661 << "SolutionHistory\n"
662 << " or getNominalValues()!\n");
663
664 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
665 initialState->setX(x);
666 }
667
668 // Perform IC Consistency
669 std::string icConsistency = this->getICConsistency();
670 TEUCHOS_TEST_FOR_EXCEPTION(
671 icConsistency != "None", std::logic_error,
672 "Error - setInitialConditions() requested a consistency of '"
673 << icConsistency
674 << "'.\n But only 'None' is available for IMEX-RK!\n");
675
676 TEUCHOS_TEST_FOR_EXCEPTION(
677 this->getUseFSAL(), std::logic_error,
678 "Error - The First-Same-As-Last (FSAL) principle is not "
679 << "available for IMEX-RK. Set useFSAL=false.\n");
680}
681
682template <typename Scalar>
684 const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& X, Scalar time,
685 Scalar stepSize, Scalar stageNumber,
686 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& G) const
687{
688 typedef Thyra::ModelEvaluatorBase MEB;
689 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
690 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
691 this->wrapperModel_);
692 MEB::InArgs<Scalar> inArgs = wrapperModelPairIMEX->createInArgs();
693 inArgs.set_x(X);
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);
698
699 // For model evaluators whose state function f(x, x_dot, t) describes
700 // an implicit ODE, and which accept an optional x_dot input argument,
701 // make sure the latter is set to null in order to request the evaluation
702 // of a state function corresponding to the explicit ODE formulation
703 // x_dot = f(x, t)
704 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
705
706 MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->createOutArgs();
707 outArgs.set_f(G);
708
709 wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs, outArgs);
710 Thyra::Vt_S(G.ptr(), -1.0);
711}
712
713template <typename Scalar>
715 const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& X, Scalar time,
716 Scalar stepSize, Scalar stageNumber,
717 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& F) const
718{
719 typedef Thyra::ModelEvaluatorBase MEB;
720
721 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
722 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
723 this->wrapperModel_);
724 MEB::InArgs<Scalar> inArgs =
725 wrapperModelPairIMEX->getExplicitModel()->createInArgs();
726 inArgs.set_x(X);
727 if (inArgs.supports(MEB::IN_ARG_t)) inArgs.set_t(time);
728 if (inArgs.supports(MEB::IN_ARG_step_size)) inArgs.set_step_size(stepSize);
729 if (inArgs.supports(MEB::IN_ARG_stage_number))
730 inArgs.set_stage_number(stageNumber);
731
732 // For model evaluators whose state function f(x, x_dot, t) describes
733 // an implicit ODE, and which accept an optional x_dot input argument,
734 // make sure the latter is set to null in order to request the evaluation
735 // of a state function corresponding to the explicit ODE formulation
736 // x_dot = f(x, t)
737 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
738
739 MEB::OutArgs<Scalar> outArgs =
740 wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
741 outArgs.set_f(F);
742
743 wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
744 Thyra::Vt_S(F.ptr(), -1.0);
745}
746
747template <class Scalar>
749 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
750{
751 this->checkInitialized();
752
753 using Teuchos::RCP;
754 using Teuchos::SerialDenseMatrix;
755 using Teuchos::SerialDenseVector;
756
757 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperIMEX_RK::takeStep()");
758 {
759 TEUCHOS_TEST_FOR_EXCEPTION(
760 solutionHistory->getNumStates() < 2, std::logic_error,
761 "Error - StepperIMEX_RK<Scalar>::takeStep(...)\n"
762 << "Need at least two SolutionStates for IMEX_RK.\n"
763 << " Number of States = " << solutionHistory->getNumStates()
764 << "\nTry setting in \"Solution History\" "
765 << "\"Storage Type\" = \"Undo\"\n"
766 << " or \"Storage Type\" = \"Static\" and "
767 << "\"Storage Limit\" = \"2\"\n");
768
769 RCP<SolutionState<Scalar> > currentState =
770 solutionHistory->getCurrentState();
771 RCP<SolutionState<Scalar> > workingState =
772 solutionHistory->getWorkingState();
773 const Scalar dt = workingState->getTimeStep();
774 const Scalar time = currentState->getTime();
775
776 const int numStages = explicitTableau_->numStages();
777 const SerialDenseMatrix<int, Scalar>& AHat = explicitTableau_->A();
778 const SerialDenseVector<int, Scalar>& bHat = explicitTableau_->b();
779 const SerialDenseVector<int, Scalar>& cHat = explicitTableau_->c();
780 const SerialDenseMatrix<int, Scalar>& A = implicitTableau_->A();
781 const SerialDenseVector<int, Scalar>& b = implicitTableau_->b();
782 const SerialDenseVector<int, Scalar>& c = implicitTableau_->c();
783
784 bool pass = true;
785 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
786
787 RCP<StepperIMEX_RK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
788 this->stepperRKAppAction_->execute(
789 solutionHistory, thisStepper,
791
792 // Compute stage solutions
793 for (int i = 0; i < numStages; ++i) {
794 this->setStageNumber(i);
795 Thyra::assign(xTilde_.ptr(), *(currentState->getX()));
796 for (int j = 0; j < i; ++j) {
797 if (AHat(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
798 Thyra::Vp_StV(xTilde_.ptr(), -dt * AHat(i, j), *(stageF_[j]));
799 if (A(i, j) != Teuchos::ScalarTraits<Scalar>::zero())
800 Thyra::Vp_StV(xTilde_.ptr(), -dt * A(i, j), *(stageG_[j]));
801 }
802
803 this->stepperRKAppAction_->execute(
804 solutionHistory, thisStepper,
806
807 Scalar ts = time + c(i) * dt;
808 Scalar tHats = time + cHat(i) * dt;
809 if (A(i, i) == Teuchos::ScalarTraits<Scalar>::zero()) {
810 // Explicit stage for the ImplicitODE_DAE
811 bool isNeeded = false;
812 for (int k = i + 1; k < numStages; ++k)
813 if (A(k, i) != 0.0) isNeeded = true;
814 if (b(i) != 0.0) isNeeded = true;
815 if (isNeeded == false) {
816 // stageG_[i] is not needed.
817 assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
818 }
819 else {
820 Thyra::assign(workingState->getX().ptr(), *xTilde_);
821 evalImplicitModelExplicitly(workingState->getX(), ts, dt, i,
822 stageG_[i]);
823 }
824 }
825 else {
826 // Implicit stage for the ImplicitODE_DAE
827 const Scalar alpha = Scalar(1.0) / (dt * A(i, i));
828 const Scalar beta = Scalar(1.0);
829
830 // Setup TimeDerivative
831 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
833 alpha, xTilde_.getConst()));
834
835 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
836 timeDer, dt, alpha, beta, SOLVE_FOR_X, i));
837
838 this->stepperRKAppAction_->execute(
839 solutionHistory, thisStepper,
841
842 const Thyra::SolveStatus<Scalar> sStatus =
843 this->solveImplicitODE(workingState->getX(), stageG_[i], ts, p);
844
845 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass = false;
846
847 this->stepperRKAppAction_->execute(
848 solutionHistory, thisStepper,
850
851 // Update contributions to stage values
852 Thyra::V_StVpStV(stageG_[i].ptr(), -alpha, *workingState->getX(), alpha,
853 *xTilde_);
854 }
855
856 this->stepperRKAppAction_->execute(
857 solutionHistory, thisStepper,
859 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
860 this->stepperRKAppAction_->execute(
861 solutionHistory, thisStepper,
863 }
864
865 // reset the stage number
866 this->setStageNumber(-1);
867
868 // Sum for solution: x_n = x_n-1 - dt*Sum{ bHat(i)*f(i) + b(i)*g(i) }
869 Thyra::assign((workingState->getX()).ptr(), *(currentState->getX()));
870 for (int i = 0; i < numStages; ++i) {
871 if (bHat(i) != Teuchos::ScalarTraits<Scalar>::zero())
872 Thyra::Vp_StV((workingState->getX()).ptr(), -dt * bHat(i),
873 *(stageF_[i]));
874 if (b(i) != Teuchos::ScalarTraits<Scalar>::zero())
875 Thyra::Vp_StV((workingState->getX()).ptr(), -dt * b(i), *(stageG_[i]));
876 }
877
878 if (pass == true)
879 workingState->setSolutionStatus(Status::PASSED);
880 else
881 workingState->setSolutionStatus(Status::FAILED);
882 workingState->setOrder(this->getOrder());
883 workingState->computeNorms(currentState);
884 this->stepperRKAppAction_->execute(
885 solutionHistory, thisStepper,
887 }
888 return;
889}
890
897template <class Scalar>
898Teuchos::RCP<Tempus::StepperState<Scalar> >
900{
901 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
902 rcp(new StepperState<Scalar>(this->getStepperType()));
903 return stepperState;
904}
905
906template <class Scalar>
908 Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const
909{
910 out.setOutputToRootOnly(0);
911
912 out << std::endl;
913 Stepper<Scalar>::describe(out, verbLevel);
914 StepperImplicit<Scalar>::describe(out, verbLevel);
915
916 out << "--- StepperIMEX_RK_Partition ---\n";
917 out << " explicitTableau_ = " << explicitTableau_ << std::endl;
918 if (verbLevel == Teuchos::VERB_HIGH)
919 explicitTableau_->describe(out, verbLevel);
920 out << " implicitTableau_ = " << implicitTableau_ << std::endl;
921 if (verbLevel == Teuchos::VERB_HIGH)
922 implicitTableau_->describe(out, verbLevel);
923 out << " xTilde_ = " << xTilde_ << std::endl;
924 out << " stageF_.size() = " << stageF_.size() << std::endl;
925 int numStages = stageF_.size();
926 for (int i = 0; i < numStages; ++i)
927 out << " stageF_[" << i << "] = " << stageF_[i] << std::endl;
928 out << " stageG_.size() = " << stageG_.size() << std::endl;
929 numStages = stageG_.size();
930 for (int i = 0; i < numStages; ++i)
931 out << " stageG_[" << i << "] = " << stageG_[i] << std::endl;
932 out << " stepperRKAppAction_= " << this->stepperRKAppAction_ << std::endl;
933 out << " order_ = " << order_ << std::endl;
934 out << "--------------------------------" << std::endl;
935}
936
937template <class Scalar>
938bool StepperIMEX_RK<Scalar>::isValidSetup(Teuchos::FancyOStream& out) const
939{
940 out.setOutputToRootOnly(0);
941 bool isValidSetup = true;
942
943 if (!Stepper<Scalar>::isValidSetup(out)) isValidSetup = false;
944
945 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
946 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
947 this->wrapperModel_);
948
949 if (wrapperModelPairIMEX->getExplicitModel() == Teuchos::null) {
950 isValidSetup = false;
951 out << "The explicit ModelEvaluator is not set!\n";
952 }
953
954 if (wrapperModelPairIMEX->getImplicitModel() == Teuchos::null) {
955 isValidSetup = false;
956 out << "The implicit ModelEvaluator is not set!\n";
957 }
958
959 if (this->wrapperModel_ == Teuchos::null) {
960 isValidSetup = false;
961 out << "The wrapper ModelEvaluator is not set!\n";
962 }
963
964 if (this->stepperRKAppAction_ == Teuchos::null) {
965 isValidSetup = false;
966 out << "The AppAction is not set!\n";
967 }
968
969 if (explicitTableau_ == Teuchos::null) {
970 isValidSetup = false;
971 out << "The explicit tableau is not set!\n";
972 }
973
974 if (implicitTableau_ == Teuchos::null) {
975 isValidSetup = false;
976 out << "The implicit tableau is not set!\n";
977 }
978
979 return isValidSetup;
980}
981
982template <class Scalar>
983Teuchos::RCP<const Teuchos::ParameterList>
985{
986 auto pl = this->getValidParametersBasicImplicit();
987 pl->template set<int>("overall order", this->getOrder());
988
989 auto explicitStepper = Teuchos::rcp(new StepperERK_General<Scalar>());
990 explicitStepper->setTableau(
991 explicitTableau_->A(), explicitTableau_->b(), explicitTableau_->c(),
992 explicitTableau_->order(), explicitTableau_->orderMin(),
993 explicitTableau_->orderMax(), explicitTableau_->bstar());
994 pl->set("IMEX-RK Explicit Stepper", *explicitStepper->getValidParameters());
995
996 auto implicitStepper = Teuchos::rcp(new StepperERK_General<Scalar>());
997 implicitStepper->setTableau(
998 implicitTableau_->A(), implicitTableau_->b(), implicitTableau_->c(),
999 implicitTableau_->order(), implicitTableau_->orderMin(),
1000 implicitTableau_->orderMax(), implicitTableau_->bstar());
1001 pl->set("IMEX-RK Implicit Stepper", *implicitStepper->getValidParameters());
1002
1003 return pl;
1004}
1005
1006// Nonmember constructor - ModelEvaluator and ParameterList
1007// ------------------------------------------------------------------------
1008template <class Scalar>
1009Teuchos::RCP<StepperIMEX_RK<Scalar> > createStepperIMEX_RK(
1010 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
1011 std::string stepperType, Teuchos::RCP<Teuchos::ParameterList> pl)
1012{
1013 auto stepper = Teuchos::rcp(new StepperIMEX_RK<Scalar>(stepperType));
1014 stepper->setStepperImplicitValues(pl);
1015 stepper->setTableaus(pl, stepperType);
1016
1017 if (model != Teuchos::null) {
1018 stepper->setModel(model);
1019 stepper->initialize();
1020 }
1021
1022 return stepper;
1023}
1024
1025} // namespace Tempus
1026#endif // Tempus_StepperIMEX_RK_impl_hpp
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
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.