Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
12 #include "Tempus_StepperRKBase.hpp"
15 
16 #include "../TestModels/DahlquistTestModel.hpp"
17 
18 
19 namespace Tempus_Unit_Test {
20 
21 using Teuchos::RCP;
22 using Teuchos::rcp;
23 using Teuchos::rcp_const_cast;
24 using Teuchos::rcp_dynamic_cast;
25 
26 
117 class StepperRKModifierBogackiShampineTest
118  : virtual public Tempus::StepperRKModifierBase<double>
119 {
120 public:
121 
123  StepperRKModifierBogackiShampineTest(Teuchos::FancyOStream &Out, bool &Success)
124  : out(Out), success(Success)
125  {}
126 
128  virtual ~StepperRKModifierBogackiShampineTest(){}
129 
131  virtual void modify(
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 
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 ||
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);
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);
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);
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);
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 
239 private:
240 
241  Teuchos::FancyOStream & out;
242  bool & success;
243 };
244 
245 
246 // ************************************************************
247 // ************************************************************
248 TEUCHOS_UNIT_TEST(ERK, Modifier)
249 {
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 // ************************************************************
279 TEUCHOS_UNIT_TEST(ERK, Modifier_FSAL)
280 {
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 // ************************************************************
310 TEUCHOS_UNIT_TEST(ERK, Modifier_FSAL_Failure)
311 {
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 // ************************************************************
342 TEUCHOS_UNIT_TEST(ERK, Modifier_Embedded)
343 {
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 
380 class StepperRKModifierXBogackiShampineTest
381  : virtual public Tempus::StepperRKModifierXBase<double>
382 {
383 public:
384 
386  StepperRKModifierXBogackiShampineTest(Teuchos::FancyOStream &Out, bool &Success)
387  : out(Out), success(Success)
388  {}
389 
391  virtual ~StepperRKModifierXBogackiShampineTest(){}
392 
394  virtual void modify(
396  const double time, const double dt, const int stageNumber,
398  {
399  const double relTol = 1.0e-14;
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 ||
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 
443 private:
444 
445  Teuchos::FancyOStream & out;
446  bool & success;
447 };
448 
449 
450 // ************************************************************
451 // ************************************************************
452 TEUCHOS_UNIT_TEST(ERK, ModifierX)
453 {
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 // ************************************************************
483 TEUCHOS_UNIT_TEST(ERK, ModifierX_FSAL)
484 {
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 // ************************************************************
514 TEUCHOS_UNIT_TEST(ERK, ModifierX_FSAL_Failure)
515 {
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_Unit_Test
The classic Dahlquist Test Problem.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Explicit RK Bogacki-Shampine Butcher Tableau.
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
MODIFIER_TYPE
Indicates the location of application action (see algorithm).
ACTION_LOCATION
Indicates the location of application action (see algorithm).
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
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.