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: 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_impl_hpp
10#define Tempus_StepperIMEX_RK_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 != "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"
39 " 'SSP1_111'\n"
40 " 'IMEX RK SSP2'\n"
41 " 'IMEX RK SSP3'\n"
42 " 'SSP3_332'\n"
43 " 'SSP2_222'\n"
44 " 'SSP2_222_L'\n"
45 " 'SSP2_222_A'\n"
46 " 'IMEX RK ARS 233'\n"
47 " 'ARS 233'\n"
48 " 'General IMEX RK'\n");
49
50 this->setStepperName( stepperType);
51 this->setStepperType( stepperType);
52 this->setUseFSAL( false);
53 this->setICConsistency( "None");
54 this->setICConsistencyCheck( false);
55 this->setZeroInitialGuess( false);
56
57 this->setStageNumber(-1);
58
59 this->setTableaus(stepperType);
60 this->setAppAction(Teuchos::null);
61 this->setDefaultSolver();
62}
63
64
65template<class Scalar>
67 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel,
68 const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
69 bool useFSAL,
70 std::string ICConsistency,
71 bool ICConsistencyCheck,
72 bool zeroInitialGuess,
73 const Teuchos::RCP<StepperRKAppAction<Scalar> >& stepperRKAppAction,
74 std::string stepperType,
75 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
76 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau,
77 Scalar order)
78{
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"
94 " 'SSP1_111'\n"
95 " 'IMEX RK SSP2'\n"
96 " 'IMEX RK SSP3'\n"
97 " 'SSP3_332'\n"
98 " 'SSP2_222'\n"
99 " 'SSP2_222_L'\n"
100 " 'SSP2_222_A'\n"
101 " 'IMEX RK ARS 233'\n"
102 " 'ARS 233'\n"
103 " 'General IMEX RK'\n");
104
105 this->setStepperName( stepperType);
106 this->setStepperType( stepperType);
107 this->setUseFSAL( useFSAL);
108 this->setICConsistency( ICConsistency);
109 this->setICConsistencyCheck( ICConsistencyCheck);
110 this->setZeroInitialGuess( zeroInitialGuess);
111
112 this->setStageNumber(-1);
113
114 if ( stepperType == "General IMEX RK" ) {
115 this->setExplicitTableau(explicitTableau);
116 this->setImplicitTableau(implicitTableau);
117 } else {
118 this->setTableaus(stepperType);
119 }
120 this->setOrder(order);
121 this->setAppAction(stepperRKAppAction);
122 this->setSolver(solver);
123
124 if (appModel != Teuchos::null) {
125 this->setModel(appModel);
126 this->initialize();
127 }
128}
129
130
131template<class Scalar>
132void StepperIMEX_RK<Scalar>::setTableaus(std::string stepperType,
133 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau,
134 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
135{
136 if (stepperType == "") stepperType = "IMEX RK SSP2";
137
138 if (stepperType == "IMEX RK 1st order" ) {
139 {
140 // Explicit Tableau
141 typedef Teuchos::ScalarTraits<Scalar> ST;
142 int NumStages = 2;
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();
148
149 // Fill A:
150 A(0,0) = zero; A(0,1) = zero;
151 A(1,0) = one; A(1,1) = zero;
152
153 // Fill b:
154 b(0) = one; b(1) = zero;
155
156 // Fill c:
157 c(0) = zero; c(1) = one;
158
159 int order = 1;
160
161 auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
162 "Explicit Tableau - IMEX RK 1st order",
163 A,b,c,order,order,order));
164 expTableau->setTVD(true);
165 expTableau->setTVDCoeff(2.0);
166
167 this->setExplicitTableau(expTableau);
168 }
169 {
170 // Implicit Tableau
171 typedef Teuchos::ScalarTraits<Scalar> ST;
172 int NumStages = 2;
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();
179
180 // Fill A:
181 A(0,0) = zero; A(0,1) = zero;
182 A(1,0) = zero; A(1,1) = one;
183
184 // Fill b:
185 b(0) = zero; b(1) = one;
186
187 // Fill c:
188 c(0) = zero; c(1) = one;
189
190 int order = 1;
191
192 auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
193 "Implicit Tableau - IMEX RK 1st order",
194 A,b,c,order,order,order));
195 impTableau->setTVD(true);
196 impTableau->setTVDCoeff(sspcoef);
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 {
207 // Explicit Tableau
208 typedef Teuchos::ScalarTraits<Scalar> ST;
209 const int NumStages = 1;
210 const int order = 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();
216
217 // Fill A:
218 A(0,0) = zero;
219
220 // Fill b:
221 b(0) = one;
222
223 // Fill c:
224 c(0) = zero;
225
226 auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
227 "Explicit Tableau - SSP1_111",A,b,c,order,order,order));
228 expTableau->setTVD(true);
229 expTableau->setTVDCoeff(1.0);
230
231 this->setExplicitTableau(expTableau);
232 }
233 {
234 // Implicit Tableau
235 typedef Teuchos::ScalarTraits<Scalar> ST;
236 const int NumStages = 1;
237 const int order = 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();
243
244 // Fill A:
245 A(0,0) = one;
246
247 // Fill b:
248 b(0) = one;
249
250 // Fill c:
251 c(0) = one;
252
253 auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
254 "Implicit Tableau - SSP1_111",
255 A,b,c,order,order,order));
256 impTableau->setTVD(true);
257 impTableau->setTVDCoeff(sspcoef);
258
259 this->setImplicitTableau(impTableau);
260 }
261 this->setStepperName("SSP1_111");
262 this->setStepperType("SSP1_111");
263 this->setOrder(1);
264
265 } else if (stepperType == "IMEX RK SSP2" || stepperType == "SSP2_222_L" ) {
266 // Explicit Tableau
267 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
268 this->setExplicitTableau(stepperERK->getTableau());
269
270 // Implicit Tableau
271 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
272 stepperSDIRK->setGammaType("2nd Order L-stable");
273 this->setImplicitTableau(stepperSDIRK->getTableau());
274
275 this->setStepperName("IMEX RK SSP2");
276 this->setStepperType("IMEX RK SSP2");
277 this->setOrder(2);
278 } else if (stepperType == "SSP2_222" || stepperType == "SSP2_222_A" ) {
279 // Explicit Tableau
280 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
281 this->setExplicitTableau(stepperERK->getTableau());
282
283 // Implicit Tableau
284 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
285 stepperSDIRK->setGammaType("gamma");
286 stepperSDIRK->setGamma( 0.5 );
287 this->setImplicitTableau(stepperSDIRK->getTableau());
288
289 this->setStepperName("SSP2_222");
290 this->setStepperType("SSP2_222");
291 this->setOrder(2);
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 } else if (stepperType == "IMEX RK ARS 233" || stepperType == "ARS 233" ) {
305 typedef Teuchos::ScalarTraits<Scalar> ST;
306 int NumStages = 3;
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);
314 {
315 // Explicit Tableau
316 // Fill A:
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;
320
321 // Fill b:
322 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
323
324 // Fill c:
325 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
326
327 int order = 2;
328
329 auto expTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
330 "Partition IMEX-RK Explicit Stepper",A,b,c,order,order,order));
331
332 this->setExplicitTableau(expTableau);
333 }
334 {
335 // Implicit Tableau
336 // Fill A:
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;
340
341 // Fill b:
342 b(0) = zero; b(1) = onehalf; b(2) = onehalf;
343
344 // Fill c:
345 c(0) = zero; c(1) = gamma; c(2) = one-gamma;
346
347 int order = 3;
348
349 auto impTableau = Teuchos::rcp(new RKButcherTableau<Scalar>(
350 "Partition IMEX-RK Implicit Stepper",A,b,c,order,order,order));
351
352 this->setImplicitTableau(impTableau);
353 }
354 this->setStepperName("IMEX RK ARS 233");
355 this->setStepperType("IMEX RK ARS 233");
356 this->setOrder(3);
357
358 } else if (stepperType == "General IMEX RK") {
359
360 if (explicitTableau == Teuchos::null) {
361 // Default Explicit Tableau (i.e., IMEX RK SSP2)
362 auto stepperERK = Teuchos::rcp(new StepperERK_Trapezoidal<Scalar>());
363 this->setExplicitTableau(stepperERK->getTableau());
364 } else {
365 this->setExplicitTableau(explicitTableau);
366 }
367
368 if (explicitTableau == Teuchos::null) {
369 // Default Implicit Tableau (i.e., IMEX RK SSP2)
370 auto stepperSDIRK = Teuchos::rcp(new StepperSDIRK_2Stage3rdOrder<Scalar>());
371 stepperSDIRK->setGammaType("2nd Order L-stable");
372 this->setImplicitTableau(stepperSDIRK->getTableau());
373 } else {
374 this->setImplicitTableau(implicitTableau);
375 }
376
377 this->setStepperName("General IMEX RK");
378 this->setStepperType("General IMEX RK");
379 this->setOrder(1);
380
381 } else {
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"
387 << " 'SSP1_111\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");
393 }
394
395 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau_==Teuchos::null,
396 std::runtime_error,
397 "Error - StepperIMEX_RK - Explicit tableau is null!");
398 TEUCHOS_TEST_FOR_EXCEPTION(implicitTableau_==Teuchos::null,
399 std::runtime_error,
400 "Error - StepperIMEX_RK - Implicit tableau is null!");
401 TEUCHOS_TEST_FOR_EXCEPTION(
402 explicitTableau_->numStages()!=implicitTableau_->numStages(),
403 std::runtime_error,
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");
409
410 this->isInitialized_ = false;
411}
412
413
414template<class Scalar>
415void
417setTableaus(Teuchos::RCP<Teuchos::ParameterList> pl,
418 std::string stepperType)
419{
420 using Teuchos::RCP;
421 if (stepperType == "") {
422 if (pl == Teuchos::null)
423 stepperType = "IMEX RK SSP2";
424 else
425 stepperType = pl->get<std::string>("Stepper Type", "IMEX RK SSP2");
426 }
427
428 if (stepperType != "General IMEX RK") {
429 this->setTableaus(stepperType);
430 } else {
431 if (pl != Teuchos::null) {
432 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau;
433 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau;
434 if (pl->isSublist("IMEX-RK Explicit Stepper")) {
435 RCP<Teuchos::ParameterList> explicitPL = Teuchos::rcp(
436 new Teuchos::ParameterList(pl->sublist("IMEX-RK Explicit Stepper")));
437 auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
438 auto stepperTemp = sf->createStepper(explicitPL, Teuchos::null);
439 auto stepperERK = Teuchos::rcp_dynamic_cast<StepperExplicitRK<Scalar> > (
440 stepperTemp,true);
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();
444 }
445
446 if (pl->isSublist("IMEX-RK Implicit Stepper")) {
447 RCP<Teuchos::ParameterList> implicitPL = Teuchos::rcp(
448 new Teuchos::ParameterList(pl->sublist("IMEX-RK Implicit Stepper")));
449 auto sf = Teuchos::rcp(new StepperFactory<Scalar>());
450 auto stepperTemp = sf->createStepper(implicitPL, Teuchos::null);
451 auto stepperDIRK = Teuchos::rcp_dynamic_cast<StepperDIRK<Scalar> > (
452 stepperTemp,true);
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();
456 }
457
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");
462
463 this->setTableaus(stepperType, explicitTableau, implicitTableau);
464 this->setOrder(pl->get<int>("overall order", 1));
465 }
466 }
467}
468
469
470template<class Scalar>
472 Teuchos::RCP<const RKButcherTableau<Scalar> > explicitTableau)
473{
474 TEUCHOS_TEST_FOR_EXCEPTION(explicitTableau->isImplicit() == true,
475 std::logic_error,
476 "Error - Received an implicit Tableau for setExplicitTableau()!\n" <<
477 " Tableau = " << explicitTableau->description() << "\n");
478 explicitTableau_ = explicitTableau;
479
480 this->isInitialized_ = false;
481}
482
483
484template<class Scalar>
486 Teuchos::RCP<const RKButcherTableau<Scalar> > implicitTableau)
487{
488 TEUCHOS_TEST_FOR_EXCEPTION( implicitTableau->isDIRK() != true,
489 std::logic_error,
490 "Error - Did not receive a DIRK Tableau for setImplicitTableau()!\n" <<
491 " Tableau = " << implicitTableau->description() << "\n");
492 implicitTableau_ = implicitTableau;
493
494 this->isInitialized_ = false;
495}
496
497template<class Scalar>
499 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& appModel)
500{
501 using Teuchos::RCP;
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");
514
515 setModelPair(modelPairIMEX);
516
517 TEUCHOS_TEST_FOR_EXCEPTION(this->solver_ == Teuchos::null, std::logic_error,
518 "Error - Solver is not set!\n");
519 this->solver_->setModel(modelPairIMEX);
520
521 this->isInitialized_ = false;
522}
523
524
530template<class Scalar>
532 const Teuchos::RCP<WrapperModelEvaluatorPairIMEX_Basic<Scalar> > &
533 modelPairIMEX)
534{
535 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
536 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
537 this->wrapperModel_);
538 validExplicitODE (modelPairIMEX->getExplicitModel());
539 validImplicitODE_DAE(modelPairIMEX->getImplicitModel());
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,
545 "Error - \n"
546 " Explicit and Implicit x vectors are incompatible!\n"
547 " Explicit vector dim = " << expXDim << "\n"
548 " Implicit vector dim = " << impXDim << "\n");
549
550 this->wrapperModel_ = wrapperModelPairIMEX;
551
552 this->isInitialized_ = false;
553}
554
560template<class Scalar>
562 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& explicitModel,
563 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& implicitModel)
564{
565 validExplicitODE (explicitModel);
566 validImplicitODE_DAE(implicitModel);
567 this->wrapperModel_ = Teuchos::rcp(
568 new WrapperModelEvaluatorPairIMEX_Basic<Scalar>(
569 explicitModel, implicitModel));
570
571 this->isInitialized_ = false;
572}
573
574
575template<class Scalar>
577{
578 TEUCHOS_TEST_FOR_EXCEPTION( this->wrapperModel_ == Teuchos::null,
579 std::logic_error,
580 "Error - Need to set the model, setModel(), before calling "
581 "StepperIMEX_RK::initialize()\n");
582
583 // Initialize the stage vectors
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());
592 }
593
594 xTilde_ = Thyra::createMember(this->wrapperModel_->get_x_space());
595 assign(xTilde_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
596
598}
599
600
601template<class Scalar>
603 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
604{
605 using Teuchos::RCP;
606
607 int numStates = solutionHistory->getNumStates();
608
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);
612
613 if (numStates > 1) {
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;
618 }
619
620 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
621 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
622
623 // Use x from inArgs as ICs, if needed.
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");
630
631 x = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
632 initialState->setX(x);
633 }
634
635 // Perform IC Consistency
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");
641
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");
645}
646
647
648template <typename Scalar>
650 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & X,
651 Scalar time, Scalar stepSize, Scalar stageNumber,
652 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & G) const
653{
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();
659 inArgs.set_x(X);
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);
664
665 // For model evaluators whose state function f(x, x_dot, t) describes
666 // an implicit ODE, and which accept an optional x_dot input argument,
667 // make sure the latter is set to null in order to request the evaluation
668 // of a state function corresponding to the explicit ODE formulation
669 // x_dot = f(x, t)
670 if (inArgs.supports(MEB::IN_ARG_x_dot)) inArgs.set_x_dot(Teuchos::null);
671
672 MEB::OutArgs<Scalar> outArgs = wrapperModelPairIMEX->getOutArgs();
673 outArgs.set_f(G);
674
675 wrapperModelPairIMEX->getImplicitModel()->evalModel(inArgs,outArgs);
676 Thyra::Vt_S(G.ptr(), -1.0);
677}
678
679
680template <typename Scalar>
682 const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & X,
683 Scalar time, Scalar stepSize, Scalar stageNumber,
684 const Teuchos::RCP<Thyra::VectorBase<Scalar> > & F) const
685{
686 typedef Thyra::ModelEvaluatorBase MEB;
687
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();
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 =
707 wrapperModelPairIMEX->getExplicitModel()->createOutArgs();
708 outArgs.set_f(F);
709
710 wrapperModelPairIMEX->getExplicitModel()->evalModel(inArgs, outArgs);
711 Thyra::Vt_S(F.ptr(), -1.0);
712}
713
714
715template<class Scalar>
717 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
718{
719 this->checkInitialized();
720
721 using Teuchos::RCP;
722 using Teuchos::SerialDenseMatrix;
723 using Teuchos::SerialDenseVector;
724
725 TEMPUS_FUNC_TIME_MONITOR("Tempus::StepperIMEX_RK::takeStep()");
726 {
727 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
728 std::logic_error,
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");
734
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();
739
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();
747
748 bool pass = true;
749 Thyra::assign(workingState->getX().ptr(), *(currentState->getX()));
750
751 RCP<StepperIMEX_RK<Scalar> > thisStepper = Teuchos::rcpFromRef(*this);
752 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
754
755 // Compute stage solutions
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]));
764 }
765
766 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
768
769 Scalar ts = time + c(i)*dt;
770 Scalar tHats = time + cHat(i)*dt;
771 if (A(i,i) == Teuchos::ScalarTraits<Scalar>::zero()) {
772 // Explicit stage for the ImplicitODE_DAE
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) {
777 // stageG_[i] is not needed.
778 assign(stageG_[i].ptr(), Teuchos::ScalarTraits<Scalar>::zero());
779 } else {
780 Thyra::assign(workingState->getX().ptr(), *xTilde_);
781 evalImplicitModelExplicitly(workingState->getX(), ts, dt, i, stageG_[i]);
782 }
783 } else {
784 // Implicit stage for the ImplicitODE_DAE
785 const Scalar alpha = Scalar(1.0)/(dt*A(i,i));
786 const Scalar beta = Scalar(1.0);
787
788 // Setup TimeDerivative
789 Teuchos::RCP<TimeDerivative<Scalar> > timeDer =
790 Teuchos::rcp(new StepperIMEX_RKTimeDerivative<Scalar>(
791 alpha, xTilde_.getConst()));
792
793 auto p = Teuchos::rcp(new ImplicitODEParameters<Scalar>(
794 timeDer, dt, alpha, beta, SOLVE_FOR_X, i));
795
796 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
798
799 const Thyra::SolveStatus<Scalar> sStatus =
800 this->solveImplicitODE(workingState->getX(), stageG_[i], ts, p);
801
802 if (sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) pass = false;
803
804 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
806
807 // Update contributions to stage values
808 Thyra::V_StVpStV(stageG_[i].ptr(), -alpha, *workingState->getX(), alpha, *xTilde_);
809 }
810
811 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
813 evalExplicitModel(workingState->getX(), tHats, dt, i, stageF_[i]);
814 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
816 }
817
818 // reset the stage number
819 this->setStageNumber(-1);
820
821 // Sum for solution: x_n = x_n-1 - dt*Sum{ bHat(i)*f(i) + b(i)*g(i) }
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]));
828 }
829
830 if (pass == true) workingState->setSolutionStatus(Status::PASSED);
831 else workingState->setSolutionStatus(Status::FAILED);
832 workingState->setOrder(this->getOrder());
833 workingState->computeNorms(currentState);
834 this->stepperRKAppAction_->execute(solutionHistory, thisStepper,
836 }
837 return;
838}
839
846template<class Scalar>
847Teuchos::RCP<Tempus::StepperState<Scalar> >
850{
851 Teuchos::RCP<Tempus::StepperState<Scalar> > stepperState =
852 rcp(new StepperState<Scalar>(this->getStepperType()));
853 return stepperState;
854}
855
856
857template<class Scalar>
859 Teuchos::FancyOStream &out,
860 const Teuchos::EVerbosityLevel verbLevel) const
861{
862 out.setOutputToRootOnly(0);
863
864 out << std::endl;
865 Stepper<Scalar>::describe(out, verbLevel);
866 StepperImplicit<Scalar>::describe(out, verbLevel);
867
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;
887}
888
889
890template<class Scalar>
891bool StepperIMEX_RK<Scalar>::isValidSetup(Teuchos::FancyOStream & out) const
892{
893 out.setOutputToRootOnly(0);
894 bool isValidSetup = true;
895
896 if ( !Stepper<Scalar>::isValidSetup(out) ) isValidSetup = false;
897
898 Teuchos::RCP<WrapperModelEvaluatorPairIMEX<Scalar> > wrapperModelPairIMEX =
899 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorPairIMEX<Scalar> >(
900 this->wrapperModel_);
901
902 if ( wrapperModelPairIMEX->getExplicitModel() == Teuchos::null) {
903 isValidSetup = false;
904 out << "The explicit ModelEvaluator is not set!\n";
905 }
906
907 if ( wrapperModelPairIMEX->getImplicitModel() == Teuchos::null) {
908 isValidSetup = false;
909 out << "The implicit ModelEvaluator is not set!\n";
910 }
911
912 if (this->wrapperModel_ == Teuchos::null) {
913 isValidSetup = false;
914 out << "The wrapper ModelEvaluator is not set!\n";
915 }
916
917 if (this->stepperRKAppAction_ == Teuchos::null) {
918 isValidSetup = false;
919 out << "The AppAction is not set!\n";
920 }
921
922 if ( explicitTableau_ == Teuchos::null ) {
923 isValidSetup = false;
924 out << "The explicit tableau is not set!\n";
925 }
926
927 if ( implicitTableau_ == Teuchos::null ) {
928 isValidSetup = false;
929 out << "The implicit tableau is not set!\n";
930 }
931
932 return isValidSetup;
933}
934
935
936template<class Scalar>
937Teuchos::RCP<const Teuchos::ParameterList>
939{
940 auto pl = this->getValidParametersBasicImplicit();
941 pl->template set<int>("overall order", this->getOrder());
942
943 auto explicitStepper = Teuchos::rcp(new StepperERK_General<Scalar>());
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());
949
950 auto implicitStepper = Teuchos::rcp(new StepperERK_General<Scalar>());
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());
956
957 return pl;
958}
959
960
961// Nonmember constructor - ModelEvaluator and ParameterList
962// ------------------------------------------------------------------------
963template<class Scalar>
964Teuchos::RCP<StepperIMEX_RK<Scalar> >
966 const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
967 std::string stepperType,
968 Teuchos::RCP<Teuchos::ParameterList> pl)
969{
970
971 auto stepper = Teuchos::rcp(new StepperIMEX_RK<Scalar>(stepperType));
972 stepper->setStepperImplicitValues(pl);
973 stepper->setTableaus(pl, stepperType);
974
975 if (model != Teuchos::null) {
976 stepper->setModel(model);
977 stepper->initialize();
978 }
979
980 return stepper;
981}
982
983
984} // namespace Tempus
985#endif // Tempus_StepperIMEX_RK_impl_hpp
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.
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< 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.