Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_UnitTest_OperatorSplit.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 
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Thyra_VectorStdOps.hpp"
15 
16 #include "Tempus_StepperFactory.hpp"
19 
20 #include "Tempus_StepperForwardEuler.hpp"
21 #include "Tempus_StepperBackwardEuler.hpp"
22 
23 #include "Tempus_StepperOperatorSplit.hpp"
30 
31 #include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
32 #include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
33 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
34 
35 #include <fstream>
36 #include <vector>
37 
38 namespace Tempus_Unit_Test {
39 
40 using Teuchos::RCP;
41 using Teuchos::rcp;
42 using Teuchos::rcp_const_cast;
43 using Teuchos::rcp_dynamic_cast;
45 using Teuchos::sublist;
46 using Teuchos::getParametersFromXmlFile;
47 
49 
50 
51 // ************************************************************
52 // ************************************************************
53 TEUCHOS_UNIT_TEST(OperatorSplit, Default_Construction)
54 {
59 
60  // Default construction.
61  auto stepper = rcp(new Tempus::StepperOperatorSplit<double>());
62  auto subStepper1 = Tempus::createStepperForwardEuler(explicitModel, Teuchos::null);
63  auto subStepper2 = Tempus::createStepperBackwardEuler(implicitModel, Teuchos::null);
64  stepper->addStepper(subStepper1);
65  stepper->addStepper(subStepper2);
66  stepper->initialize();
67  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
68 
69  // Default values for construction.
73  bool useFSAL = stepper->getUseFSAL();
74  std::string ICConsistency = stepper->getICConsistency();
75  bool ICConsistencyCheck = stepper->getICConsistencyCheck();
76  int order = 1;
77 
78  // Test the set functions.
79  stepper->setAppAction(modifier); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
80  stepper->setAppAction(modifierX); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
81  stepper->setAppAction(observer); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
82  stepper->setUseFSAL(useFSAL); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
83  stepper->setICConsistency(ICConsistency); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
84  stepper->setICConsistencyCheck(ICConsistencyCheck); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
85  stepper->setOrder(order); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
86  stepper->setOrderMin(order); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
87  stepper->setOrderMax(order); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
88 
89 
90  // Full argument list construction.
91  std::vector<RCP<const Thyra::ModelEvaluator<double> > > models;
92  models.push_back(explicitModel);
93  models.push_back(implicitModel);
94 
95  std::vector<Teuchos::RCP<Tempus::Stepper<double> > > subStepperList;
96  subStepperList.push_back(subStepper1);
97  subStepperList.push_back(subStepper2);
98 
100  models, subStepperList, useFSAL, ICConsistency, ICConsistencyCheck,order, order, order,modifier));
101 
102  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
103 
104  // Test stepper properties.
105  TEUCHOS_ASSERT(stepper->getOrder() == 1);
106 }
107 
108 
109 // ************************************************************
110 // ************************************************************
111 TEUCHOS_UNIT_TEST(OperatorSplit, StepperFactory_Construction)
112 {
113  // Read params from .xml file
114  auto pList = getParametersFromXmlFile(
115  "../test/OperatorSplit/Tempus_OperatorSplit_VanDerPol.xml");
116  auto tempusPL = sublist(pList, "Tempus", true);
117  auto stepperPL = sublist(tempusPL, "Demo Stepper", true);
118 
119  auto explicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ExplicitModel<double>());
120  auto implicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ImplicitModel<double>());
121  std::vector<RCP<const Thyra::ModelEvaluator<double> > > models;
122  models.push_back(explicitModel);
123  models.push_back(implicitModel);
124 
125 
127 
128  // Test using ParameterList.
129  // Passing in model.
130  auto stepper = sf->createStepper(stepperPL, models);
131  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
132 
133 }
134 
135 // ************************************************************
136 // ************************************************************
137 class StepperOperatorSplitModifierTest
138  : virtual public Tempus::StepperOperatorSplitModifierBase<double>
139 {
140 public:
141 
143  StepperOperatorSplitModifierTest()
144  : testBEGIN_STEP(false), testEND_STEP(false),
145  testCurrentValue(-0.99), testWorkingValue(-0.99),
146  testDt(-1.5), testName("")
147  {}
148 
150  virtual ~StepperOperatorSplitModifierTest(){}
151 
153  virtual void modify(
157  {
158  switch(actLoc) {
160  {
161  testBEGIN_STEP = true;
162  auto x = sh->getCurrentState()->getX();
163  testCurrentValue = get_ele(*(x), 0);
164  break;
165  }
167  {
168  testBEFORE_STEPPER = true;
169  testDt = sh->getWorkingState()->getTimeStep()/10.0;
170  sh->getWorkingState()->setTimeStep(testDt);
171  break;
172  }
174  {
175  testAFTER_STEPPER = true;
176  testName = "OperatorSplit - Modifier";
177  stepper->setStepperName(testName);
178  break;
179  }
181  {
182  testEND_STEP = true;
183  auto x = sh->getWorkingState()->getX();
184  testWorkingValue = get_ele(*(x), 0);
185  break;
186  }
187  default:
188  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
189  "Error - unknown action location.\n");
190  }
191  }
192 
193  bool testBEGIN_STEP;
194  bool testBEFORE_STEPPER;
195  bool testAFTER_STEPPER;
196  bool testEND_STEP;
197  double testCurrentValue;
198  double testWorkingValue;
199  double testDt;
200  std::string testName;
201 };
202 
203 TEUCHOS_UNIT_TEST(OperatorSplit, AppAction_Modifier)
204 {
209  // Default construction.
210  auto stepper = rcp(new Tempus::StepperOperatorSplit<double>());
211  auto subStepper1 = Tempus::createStepperForwardEuler(explicitModel, Teuchos::null);
212  auto subStepper2 = Tempus::createStepperBackwardEuler(implicitModel, Teuchos::null);
213  auto modifier = rcp(new StepperOperatorSplitModifierTest());
214  stepper->setAppAction(modifier);
215  stepper->addStepper(subStepper1);
216  stepper->addStepper(subStepper2);
217  stepper->initialize();
218 
219  // Setup initial condition SolutionState --------------------
220  auto inArgsIC = stepper->getModel()->getNominalValues();
221  auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
222  auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
223  auto icState = Tempus::createSolutionStateX(icX, icXDot);
224  icState->setTime (0.0);
225  icState->setIndex (1);
226  icState->setTimeStep(-15.0);
227  icState->setOrder (stepper->getOrder());
228  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
229 
230  // Create a SolutionHistory.
231  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
232  solutionHistory->setName("Forward States");
233  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
234  solutionHistory->setStorageLimit(2);
235  solutionHistory->addState(icState);
236 
237  // Take one time step.
238  stepper->setInitialConditions(solutionHistory);
239  solutionHistory->initWorkingState();
240  solutionHistory->getWorkingState()->setTimeStep(-15.0);
241  stepper->takeStep(solutionHistory);
242 
243  // Testing that each ACTION_LOCATION has been called.
244  TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
245  TEST_COMPARE(modifier->testBEFORE_STEPPER, ==, true);
246  TEST_COMPARE(modifier->testAFTER_STEPPER, ==, true);
247  TEST_COMPARE(modifier->testEND_STEP, ==, true);
248 
249  // Testing that values can be set through the Modifier.
250  auto x = solutionHistory->getCurrentState()->getX();
251  TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
252  x = solutionHistory->getWorkingState()->getX();
253  TEST_FLOATING_EQUALITY(modifier->testWorkingValue, get_ele(*(x), 0), 1.0e-14);
254  auto Dt = solutionHistory->getWorkingState()->getTimeStep();
255  TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
256 
257  TEST_COMPARE(modifier->testName, ==, "OperatorSplit - Modifier");
258 }
259 
260 // ************************************************************
261 // ************************************************************
262 class StepperOperatorSplitObserverTest
263  : virtual public Tempus::StepperOperatorSplitObserverBase<double>
264 {
265 public:
266 
268  StepperOperatorSplitObserverTest()
269  : testBEGIN_STEP(false), testBEFORE_STEPPER(false),
270  testAFTER_STEPPER(false), testEND_STEP(false),
271  testCurrentValue(-0.99), testWorkingValue(-0.99),
272  testDt(-1.5), testName("Operator Split")
273  {}
274 
276  virtual ~StepperOperatorSplitObserverTest(){}
277 
279  virtual void observe(
283  {
284  switch(actLoc) {
286  {
287  testBEGIN_STEP = true;
288  auto x = sh->getCurrentState()->getX();
289  testCurrentValue = get_ele(*(x), 0);
290  break;
291  }
293  {
294  testBEFORE_STEPPER = true;
295  testDt = sh->getWorkingState()->getTimeStep();
296  break;
297  }
299  {
300  testAFTER_STEPPER = true;
301  testName = stepper->getStepperType();
302  break;
303  }
305  {
306  testEND_STEP = true;
307  auto x = sh->getWorkingState()->getX();
308  testWorkingValue = get_ele(*(x), 0);
309  break;
310  }
311  default:
312  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
313  "Error - unknown action location.\n");
314  }
315  }
316 
317  bool testBEGIN_STEP;
318  bool testBEFORE_STEPPER;
319  bool testAFTER_STEPPER;
320  bool testEND_STEP;
321  double testCurrentValue;
322  double testWorkingValue;
323  double testDt;
324  std::string testName;
325 };
326 
327 TEUCHOS_UNIT_TEST(OperatorSplit, AppAction_Observer)
328 {
333  // Default construction.
334  auto stepper = rcp(new Tempus::StepperOperatorSplit<double>());
335  auto subStepper1 = Tempus::createStepperForwardEuler(explicitModel, Teuchos::null);
336  auto subStepper2 = Tempus::createStepperBackwardEuler(implicitModel, Teuchos::null);
337  auto observer = rcp(new StepperOperatorSplitObserverTest());
338  stepper->setAppAction(observer);
339  stepper->addStepper(subStepper1);
340  stepper->addStepper(subStepper2);
341  stepper->initialize();
342 
343  // Setup initial condition SolutionState --------------------
344  auto inArgsIC = stepper->getModel()->getNominalValues();
345  auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
346  auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
347  auto icState = Tempus::createSolutionStateX(icX, icXDot);
348  icState->setTime (0.0);
349  icState->setIndex (1);
350  icState->setTimeStep(-1.5);
351  icState->setOrder (stepper->getOrder());
352  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
353 
354  // Create a SolutionHistory.
355  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
356  solutionHistory->setName("Forward States");
357  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
358  solutionHistory->setStorageLimit(2);
359  solutionHistory->addState(icState);
360 
361  // Take one time step.
362  stepper->setInitialConditions(solutionHistory);
363  solutionHistory->initWorkingState();
364  solutionHistory->getWorkingState()->setTimeStep(-1.5);
365  stepper->takeStep(solutionHistory);
366 
367  // Testing that each ACTION_LOCATION has been called.
368  TEST_COMPARE(observer->testBEGIN_STEP, ==, true);
369  TEST_COMPARE(observer->testBEFORE_STEPPER, ==, true);
370  TEST_COMPARE(observer->testAFTER_STEPPER, ==, true);
371  TEST_COMPARE(observer->testEND_STEP, ==, true);
372 
373  // Testing that values can be observed through the observer.
374  auto x = solutionHistory->getCurrentState()->getX();
375  TEST_FLOATING_EQUALITY(observer->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
376  x = solutionHistory->getWorkingState()->getX();
377  TEST_FLOATING_EQUALITY(observer->testWorkingValue, get_ele(*(x), 0), 1.0e-14);
378  TEST_FLOATING_EQUALITY(observer->testDt, -1.5, 1.0e-14);
379 
380  TEST_COMPARE(observer->testName, ==, "Operator Split");
381 }
382 
383 // ************************************************************
384 // ************************************************************
385 class StepperOperatorSplitModifierXTest
386  : virtual public Tempus::StepperOperatorSplitModifierXBase<double>
387 {
388 public:
389 
391  StepperOperatorSplitModifierXTest()
392  : testX_BEGIN_STEP(false), testX_BEFORE_STEPPER(false),
393  testX_AFTER_STEPPER(false), testXDOT_END_STEP(false),
394  testX(-0.99), testXDot(-0.99),
395  testDt(-1.5), testTime(-1.5)
396  {}
397 
399  virtual ~StepperOperatorSplitModifierXTest(){}
400 
402  virtual void modify(
404  const double time, const double dt,
406  {
407  switch(modType) {
409  {
410  testX_BEGIN_STEP = true;
411  testX = get_ele(*(x), 0);
412  break;
413  }
415  {
416  testX_BEFORE_STEPPER = true;
417  testDt = dt;
418  break;
419  }
421  {
422  testX_AFTER_STEPPER = true;
423  testTime = time;
424  break;
425  }
427  {
428  testXDOT_END_STEP = true;
429  testXDot = get_ele(*(x), 0);
430  break;
431  }
432  default:
433  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
434  "Error - unknown action location.\n");
435  }
436  }
437 
438  bool testX_BEGIN_STEP;
439  bool testX_BEFORE_STEPPER;
440  bool testX_AFTER_STEPPER;
441  bool testXDOT_END_STEP;
442  double testX;
443  double testXDot;
444  double testDt;
445  double testTime;
446 };
447 
448 TEUCHOS_UNIT_TEST(OperatorSplit, AppAction_ModifierX)
449 {
454  // Default construction.
455  auto stepper = rcp(new Tempus::StepperOperatorSplit<double>());
456  auto subStepper1 = Tempus::createStepperForwardEuler(explicitModel, Teuchos::null);
457  auto subStepper2 = Tempus::createStepperBackwardEuler(implicitModel, Teuchos::null);
458  auto modifierX = rcp(new StepperOperatorSplitModifierXTest());
459  stepper->setAppAction(modifierX);
460  stepper->addStepper(subStepper1);
461  stepper->addStepper(subStepper2);
462  stepper->initialize();
463 
464  // Setup initial condition SolutionState --------------------
465  auto inArgsIC = stepper->getModel()->getNominalValues();
466  auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
467  auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
468  auto icState = Tempus::createSolutionStateX(icX, icXDot);
469  icState->setTime (0.0);
470  icState->setIndex (1);
471  icState->setTimeStep(-1.5);
472  icState->setOrder (stepper->getOrder());
473  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
474 
475  // Create a SolutionHistory.
476  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
477  solutionHistory->setName("Forward States");
478  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
479  solutionHistory->setStorageLimit(2);
480  solutionHistory->addState(icState);
481 
482  // Take one time step.
483  stepper->setInitialConditions(solutionHistory);
484  solutionHistory->initWorkingState();
485  solutionHistory->getWorkingState()->setTimeStep(-1.5);
486  stepper->takeStep(solutionHistory);
487 
488  // Testing that each ACTION_LOCATION has been called.
489  TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
490  TEST_COMPARE(modifierX->testX_BEFORE_STEPPER, ==, true);
491  TEST_COMPARE(modifierX->testX_AFTER_STEPPER, ==, true);
492  TEST_COMPARE(modifierX->testXDOT_END_STEP, ==, true);
493 
494  // Testing that values can be set through the Modifier.
495  auto x = solutionHistory->getCurrentState()->getX();
496  TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
497  // Temporary memory for xDot is not guarranteed to exist outside the Stepper.
498  auto xDot = solutionHistory->getWorkingState()->getXDot();
499  if (xDot == Teuchos::null) xDot = stepper->getStepperXDot();
500 
501  TEST_FLOATING_EQUALITY(modifierX->testXDot, get_ele(*(xDot), 0),1.0e-14);
502  auto Dt = solutionHistory->getWorkingState()->getTimeStep();
503  TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
504 
505  auto time = solutionHistory->getWorkingState()->getTime();
506  TEST_FLOATING_EQUALITY(modifierX->testTime, time, 1.0e-14);
507 }
508 
509 } // namespace Tempus_Test
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
Explicit Runge-Kutta time stepper.
MODIFIER_TYPE
Indicates the location of application action (see algorithm).
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
OperatorSplit stepper loops through the Stepper list.
Teuchos::RCP< StepperForwardEuler< Scalar > > createStepperForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
virtual void modify(Teuchos::RCP< Thyra::VectorBase< double > >, const double, const double, const MODIFIER_TYPE modType)=0
Modify solution based on the MODIFIER_TYPE.
Teuchos::RCP< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)