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