Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_UnitTest_ERK.cpp
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
10
15
16#include "../TestModels/DahlquistTestModel.hpp"
17
18
19namespace Tempus_Unit_Test {
20
21using Teuchos::RCP;
22using Teuchos::rcp;
23using Teuchos::rcp_const_cast;
24using Teuchos::rcp_dynamic_cast;
25
26
117class StepperRKModifierBogackiShampineTest
118 : virtual public Tempus::StepperRKModifierBase<double>
119{
120public:
121
123 StepperRKModifierBogackiShampineTest(Teuchos::FancyOStream &Out, bool &Success)
124 : out(Out), success(Success)
125 {}
126
128 virtual ~StepperRKModifierBogackiShampineTest(){}
129
131 virtual void modify(
132 Teuchos::RCP<Tempus::SolutionHistory<double> > sh,
133 Teuchos::RCP<Tempus::StepperRKBase<double> > stepper,
135 {
136 const double relTol = 1.0e-14;
137 auto stageNumber = stepper->getStageNumber();
138 Teuchos::SerialDenseVector<int,double> c = stepper->getTableau()->c();
139
140 auto currentState = sh->getCurrentState();
141 auto workingState = sh->getWorkingState();
142 const double dt = workingState->getTimeStep();
143 double time = currentState->getTime();
144 if (stageNumber >= 0) time += c(stageNumber)*dt;
145
146 auto x = workingState->getX();
147 auto xDot = workingState->getXDot();
148 if (xDot == Teuchos::null) xDot = stepper->getStepperXDot();
149
150
151 if (actLoc == StepperRKAppAction<double>::BEGIN_STEP) {
152 {
153 auto DME = Teuchos::rcp_dynamic_cast<
154 const Tempus_Test::DahlquistTestModel<double> > (stepper->getModel());
155 TEST_FLOATING_EQUALITY(DME->getLambda(), -1.0, relTol);
156 }
157 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
158
159 const double x_0 = get_ele(*(x), 0);
160 const double xDot_0 = get_ele(*(xDot), 0);
161 TEST_FLOATING_EQUALITY(x_0, 1.0, relTol); // Should be x_0
162 TEST_FLOATING_EQUALITY(xDot_0, -1.0, relTol); // Should be xDot_0
163 TEST_ASSERT(std::abs(time) < relTol);
164 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
165 TEST_COMPARE(stageNumber, ==, -1);
166
167 } else if (actLoc == StepperRKAppAction<double>::BEGIN_STAGE ||
168 actLoc == StepperRKAppAction<double>::BEFORE_SOLVE ||
169 actLoc == StepperRKAppAction<double>::AFTER_SOLVE ||
170 actLoc == StepperRKAppAction<double>::BEFORE_EXPLICIT_EVAL ||
171 actLoc == StepperRKAppAction<double>::END_STAGE) {
172 const double X_i = get_ele(*(x), 0);
173 const double f_i = get_ele(*(xDot), 0);
174 if (stageNumber == 0) {
175 TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); // Should be X_1
176 TEST_ASSERT(std::abs(time) < relTol);
177 if (actLoc == StepperRKAppAction<double>::END_STAGE ||
178 !stepper->getUseFSAL()) {
179 TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); // Should be \bar{f}_1
180 } else {
181 TEST_ASSERT(std::abs(f_i) < relTol); // Not set yet.
182 }
183
184 } else if (stageNumber == 1) {
185 TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); // Should be X_2
186 TEST_FLOATING_EQUALITY(time, 0.5, relTol);
187 if (actLoc == StepperRKAppAction<double>::END_STAGE) {
188 TEST_FLOATING_EQUALITY(f_i, -0.5, relTol); // Should be \bar{f}_2
189 } else {
190 TEST_ASSERT(std::abs(f_i) < relTol); // Not set yet.
191 }
192
193 } else if (stageNumber == 2) {
194 TEST_FLOATING_EQUALITY(X_i, 5.0/8.0, relTol); // Should be X_3
195 TEST_FLOATING_EQUALITY(time, 0.75, relTol);
196 if (actLoc == StepperRKAppAction<double>::END_STAGE) {
197 TEST_FLOATING_EQUALITY(f_i, -5.0/8.0, relTol); // Should be \bar{f}_3
198 } else {
199 TEST_ASSERT(std::abs(f_i) < relTol); // Not set yet.
200 }
201
202 } else if (stageNumber == 3) {
203 TEST_FLOATING_EQUALITY(X_i, 1.0/3.0, relTol); // Should be X_4
204 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
205 if (actLoc == StepperRKAppAction<double>::END_STAGE) {
206 TEST_FLOATING_EQUALITY(f_i, -1.0/3.0, relTol); // Should be \bar{f}_4
207 } else if (workingState->getNConsecutiveFailures() > 0) {
208 TEST_FLOATING_EQUALITY(f_i, -1.0, relTol); // From IC and FSAL swap.
209 } else {
210 TEST_ASSERT(std::abs(f_i) < relTol); // Not set yet.
211 }
212
213 } else {
214 TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 4));
215 }
216 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
217
218 } else if (actLoc == StepperRKAppAction<double>::END_STEP) {
219 const double x_1 = get_ele(*(x), 0);
220 time = workingState->getTime();
221 TEST_FLOATING_EQUALITY(x_1, 1.0/3.0, relTol); // Should be x_1
222 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
223 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
224 TEST_COMPARE(stageNumber, ==, -1);
225
226 if (stepper->getUseEmbedded() == true) {
227 TEST_FLOATING_EQUALITY(workingState->getTolRel(), 1.0, relTol);
228 TEST_ASSERT(std::abs(workingState->getTolAbs()) < relTol);
229 // e = 0 from doxygen above.
230 TEST_ASSERT(std::abs(workingState->getErrorRel()) < relTol);
231 }
232
233 } else {
234 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
235 "Error - unknown action location.\n");
236 }
237 }
238
239private:
240
241 Teuchos::FancyOStream & out;
242 bool & success;
243};
244
245
246// ************************************************************
247// ************************************************************
248TEUCHOS_UNIT_TEST(ERK, Modifier)
249{
250 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
251 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
253 auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));
254
255 stepper->setModel(model);
256 stepper->setAppAction(modifier);
257 stepper->setICConsistency("Consistent");
258 stepper->setUseFSAL(false);
259 stepper->initialize();
260
261 // Create a SolutionHistory.
262 auto solutionHistory = Tempus::createSolutionHistoryME(model);
263
264 // Take one time step.
265 stepper->setInitialConditions(solutionHistory);
266 solutionHistory->initWorkingState();
267 double dt = 1.0;
268 solutionHistory->getWorkingState()->setTimeStep(dt);
269 solutionHistory->getWorkingState()->setTime(dt);
270 stepper->takeStep(solutionHistory); // Primary testing occurs in
271 // modifierX during takeStep().
272 // Test stepper properties.
273 TEUCHOS_ASSERT(stepper->getOrder() == 3);
274}
275
276
277// ************************************************************
278// ************************************************************
279TEUCHOS_UNIT_TEST(ERK, Modifier_FSAL)
280{
281 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
282 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
284 auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));
285
286 stepper->setModel(model);
287 stepper->setAppAction(modifier);
288 stepper->setICConsistency("Consistent");
289 stepper->setUseFSAL(true);
290 stepper->initialize();
291
292 // Create a SolutionHistory.
293 auto solutionHistory = Tempus::createSolutionHistoryME(model);
294
295 // Take one time step.
296 stepper->setInitialConditions(solutionHistory);
297 solutionHistory->initWorkingState();
298 double dt = 1.0;
299 solutionHistory->getWorkingState()->setTimeStep(dt);
300 solutionHistory->getWorkingState()->setTime(dt);
301 stepper->takeStep(solutionHistory); // Primary testing occurs in
302 // modifierX during takeStep().
303 // Test stepper properties.
304 TEUCHOS_ASSERT(stepper->getOrder() == 3);
305}
306
307
308// ************************************************************
309// ************************************************************
310TEUCHOS_UNIT_TEST(ERK, Modifier_FSAL_Failure)
311{
312 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
313 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
315 auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));
316
317 stepper->setModel(model);
318 stepper->setAppAction(modifier);
319 stepper->setICConsistency("Consistent");
320 stepper->setUseFSAL(true);
321 stepper->initialize();
322
323 // Create a SolutionHistory.
324 auto solutionHistory = Tempus::createSolutionHistoryME(model);
325
326 // Take one time step.
327 stepper->setInitialConditions(solutionHistory);
328 solutionHistory->initWorkingState();
329 double dt = 1.0;
330 solutionHistory->getWorkingState()->setTimeStep(dt);
331 solutionHistory->getWorkingState()->setTime(dt);
332 solutionHistory->getWorkingState()->setNConsecutiveFailures(1);
333 stepper->takeStep(solutionHistory); // Primary testing occurs in
334 // modifierX during takeStep().
335 // Test stepper properties.
336 TEUCHOS_ASSERT(stepper->getOrder() == 3);
337}
338
339
340// ************************************************************
341// ************************************************************
342TEUCHOS_UNIT_TEST(ERK, Modifier_Embedded)
343{
344 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
345 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
347 auto modifier = rcp(new StepperRKModifierBogackiShampineTest(out, success));
348
349 stepper->setModel(model);
350 stepper->setAppAction(modifier);
351 stepper->setICConsistency("Consistent");
352 stepper->setUseFSAL(false);
353 stepper->setUseEmbedded(true);
354 stepper->initialize();
355
356 // Create a SolutionHistory.
357 auto solutionHistory = Tempus::createSolutionHistoryME(model);
358
359 // Take one time step.
360 stepper->setInitialConditions(solutionHistory);
361 solutionHistory->initWorkingState();
362 double dt = 1.0;
363 solutionHistory->getWorkingState()->setTimeStep(dt);
364 solutionHistory->getWorkingState()->setTime(dt);
365 solutionHistory->getWorkingState()->setTolRel(1.0);
366 solutionHistory->getWorkingState()->setTolAbs(0.0);
367 stepper->takeStep(solutionHistory); // Primary testing occurs in
368 // modifierX during takeStep().
369 // Test stepper properties.
370 TEUCHOS_ASSERT(stepper->getOrder() == 3);
371}
372
373
380class StepperRKModifierXBogackiShampineTest
381 : virtual public Tempus::StepperRKModifierXBase<double>
382{
383public:
384
386 StepperRKModifierXBogackiShampineTest(Teuchos::FancyOStream &Out, bool &Success)
387 : out(Out), success(Success)
388 {}
389
391 virtual ~StepperRKModifierXBogackiShampineTest(){}
392
394 virtual void modify(
395 Teuchos::RCP<Thyra::VectorBase<double> > x,
396 const double time, const double dt, const int stageNumber,
398 {
399 const double relTol = 1.0e-14;
400 if (modType == StepperRKModifierXBase<double>::X_BEGIN_STEP) {
401 const double x_0 = get_ele(*(x), 0); // Should be x_0
402 TEST_FLOATING_EQUALITY(x_0, 1.0, relTol);
403 TEST_FLOATING_EQUALITY(time, 0.0, relTol);
404 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
405 TEST_COMPARE(stageNumber, ==, -1);
406
407 } else if (modType == StepperRKModifierXBase<double>::X_BEGIN_STAGE ||
408 modType == StepperRKModifierXBase<double>::X_BEFORE_SOLVE ||
409 modType == StepperRKModifierXBase<double>::X_AFTER_SOLVE ||
410 modType == StepperRKModifierXBase<double>::X_BEFORE_EXPLICIT_EVAL ||
411 modType == StepperRKModifierXBase<double>::X_END_STAGE) {
412 const double X_i = get_ele(*(x), 0);
413 if (stageNumber == 0) {
414 TEST_FLOATING_EQUALITY(X_i, 1.0, relTol); // Should be X_1
415 TEST_FLOATING_EQUALITY(time, 0.0, relTol);
416 } else if (stageNumber == 1) {
417 TEST_FLOATING_EQUALITY(X_i, 0.5, relTol); // Should be X_2
418 TEST_FLOATING_EQUALITY(time, 0.5, relTol);
419 } else if (stageNumber == 2) {
420 TEST_FLOATING_EQUALITY(X_i, 5.0/8.0, relTol); // Should be X_3
421 TEST_FLOATING_EQUALITY(time, 0.75, relTol);
422 } else if (stageNumber == 3) {
423 TEST_FLOATING_EQUALITY(X_i, 1.0/3.0, relTol); // Should be X_4
424 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
425 } else {
426 TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 4));
427 }
428 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
429
430 } else if (modType == StepperRKModifierXBase<double>::X_END_STEP) {
431 const double x_1 = get_ele(*(x), 0);
432 TEST_FLOATING_EQUALITY(x_1, 1.0/3.0, relTol); // Should be x_1
433 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
434 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
435 TEST_COMPARE(stageNumber, ==, -1);
436
437 } else {
438 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
439 "Error - unknown action location.\n");
440 }
441 }
442
443private:
444
445 Teuchos::FancyOStream & out;
446 bool & success;
447};
448
449
450// ************************************************************
451// ************************************************************
452TEUCHOS_UNIT_TEST(ERK, ModifierX)
453{
454 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
455 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
457 auto modifierX = rcp(new StepperRKModifierXBogackiShampineTest(out, success));
458
459 stepper->setModel(model);
460 stepper->setAppAction(modifierX);
461 stepper->setICConsistency("Consistent");
462 stepper->setUseFSAL(false);
463 stepper->initialize();
464
465 // Create a SolutionHistory.
466 auto solutionHistory = Tempus::createSolutionHistoryME(model);
467
468 // Take one time step.
469 stepper->setInitialConditions(solutionHistory);
470 solutionHistory->initWorkingState();
471 double dt = 1.0;
472 solutionHistory->getWorkingState()->setTimeStep(dt);
473 solutionHistory->getWorkingState()->setTime(dt);
474 stepper->takeStep(solutionHistory); // Primary testing occurs in
475 // modifierX during takeStep().
476 // Test stepper properties.
477 TEUCHOS_ASSERT(stepper->getOrder() == 3);
478}
479
480
481// ************************************************************
482// ************************************************************
483TEUCHOS_UNIT_TEST(ERK, ModifierX_FSAL)
484{
485 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
486 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
488 auto modifierX = rcp(new StepperRKModifierXBogackiShampineTest(out, success));
489
490 stepper->setModel(model);
491 stepper->setAppAction(modifierX);
492 stepper->setICConsistency("Consistent");
493 stepper->setUseFSAL(true);
494 stepper->initialize();
495
496 // Create a SolutionHistory.
497 auto solutionHistory = Tempus::createSolutionHistoryME(model);
498
499 // Take one time step.
500 stepper->setInitialConditions(solutionHistory);
501 solutionHistory->initWorkingState();
502 double dt = 1.0;
503 solutionHistory->getWorkingState()->setTimeStep(dt);
504 solutionHistory->getWorkingState()->setTime(dt);
505 stepper->takeStep(solutionHistory); // Primary testing occurs in
506 // modifierX during takeStep().
507 // Test stepper properties.
508 TEUCHOS_ASSERT(stepper->getOrder() == 3);
509}
510
511
512// ************************************************************
513// ************************************************************
514TEUCHOS_UNIT_TEST(ERK, ModifierX_FSAL_Failure)
515{
516 auto stepper = rcp(new Tempus::StepperERK_BogackiShampine32<double>());
517 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
519 auto modifierX = rcp(new StepperRKModifierXBogackiShampineTest(out, success));
520
521 stepper->setModel(model);
522 stepper->setAppAction(modifierX);
523 stepper->setICConsistency("Consistent");
524 stepper->setUseFSAL(true);
525 stepper->initialize();
526
527 // Create a SolutionHistory.
528 auto solutionHistory = Tempus::createSolutionHistoryME(model);
529
530 // Take one time step.
531 stepper->setInitialConditions(solutionHistory);
532 solutionHistory->initWorkingState();
533 double dt = 1.0;
534 solutionHistory->getWorkingState()->setTimeStep(dt);
535 solutionHistory->getWorkingState()->setTime(dt);
536 solutionHistory->getWorkingState()->setNConsecutiveFailures(1);
537 stepper->takeStep(solutionHistory); // Primary testing occurs in
538 // modifierX during takeStep().
539 // Test stepper properties.
540 TEUCHOS_ASSERT(stepper->getOrder() == 3);
541}
542
543
544} // namespace Tempus_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Explicit RK Bogacki-Shampine Butcher Tableau.
ACTION_LOCATION
Indicates the location of application action (see algorithm).
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
MODIFIER_TYPE
Indicates the location of application action (see algorithm).
virtual void modify(Teuchos::RCP< Thyra::VectorBase< double > >, const double, const double, const int, const MODIFIER_TYPE modType)=0
Modify solution based on the MODIFIER_TYPE.
The classic Dahlquist Test Problem.
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.