Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_UnitTest_Subcycling.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 
11 #include "Tempus_StepperForwardEuler.hpp"
12 #include "Tempus_StepperBackwardEuler.hpp"
13 
14 #include "Tempus_StepperSubcycling.hpp"
21 
22 namespace Tempus_Unit_Test {
23 
25 using Teuchos::RCP;
26 using Teuchos::rcp;
27 using Teuchos::rcp_const_cast;
28 using Teuchos::rcp_dynamic_cast;
29 using Teuchos::sublist;
30 
32 
33 // ************************************************************
34 // ************************************************************
35 TEUCHOS_UNIT_TEST(Subcycling, Default_Construction)
36 {
37  auto model = rcp(new Tempus_Test::SinCosModel<double>());
38  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double> >(model);
39 
40  // Setup SolutionHistory ------------------------------------
41  auto inArgsIC = model->getNominalValues();
42  auto icSolution =
43  rcp_const_cast<Thyra::VectorBase<double> >(inArgsIC.get_x());
44  auto icState = Tempus::createSolutionStateX(icSolution);
45  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
46  solutionHistory->addState(icState);
47  solutionHistory->initWorkingState();
48 
49  // Default construction.
50  auto stepper = rcp(new Tempus::StepperSubcycling<double>());
51  auto stepperBE = Tempus::createStepperBackwardEuler(modelME, Teuchos::null);
52  stepper->setSubcyclingStepper(stepperBE);
53  stepper->setInitialConditions(solutionHistory);
54  stepper->initialize();
55  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
56 
57  // Default values for construction.
61  auto solver = rcp(new Thyra::NOXNonlinearSolver());
62  solver->setParameterList(Tempus::defaultSolverParameters());
63 
64  bool useFSAL = stepper->getUseFSAL();
65  std::string ICConsistency = stepper->getICConsistency();
66  bool ICConsistencyCheck = stepper->getICConsistencyCheck();
67 
68  // Test the set functions
69  stepper->setSolver(solver);
70  stepper->initialize();
71  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
72  stepper->setAppAction(modifier);
73  stepper->initialize();
74  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
75  stepper->setAppAction(modifierX);
76  stepper->initialize();
77  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
78  stepper->setAppAction(observer);
79  stepper->initialize();
80  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
81  stepper->setUseFSAL(useFSAL);
82  stepper->initialize();
83  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
84  stepper->setICConsistency(ICConsistency);
85  stepper->initialize();
86  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
87  stepper->setICConsistencyCheck(ICConsistencyCheck);
88  stepper->initialize();
89  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
90 
91  // Full argument list construction.
92  auto scIntegrator = Teuchos::rcp(new Tempus::IntegratorBasic<double>());
93  auto stepperFE = Tempus::createStepperForwardEuler(modelME, Teuchos::null);
94  scIntegrator->setStepper(stepperFE);
95  scIntegrator->setSolutionHistory(solutionHistory);
96  scIntegrator->initialize();
97 
99  model, scIntegrator, useFSAL, ICConsistency, ICConsistencyCheck,
100  modifier));
101  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
102 
103  // Test stepper properties.
104  TEUCHOS_ASSERT(stepper->getOrder() == 1);
105 }
106 
107 // ************************************************************
108 // ************************************************************
109 TEUCHOS_UNIT_TEST(Subcycling, MaxTimeStepDoesNotChangeDuring_takeStep)
110 {
111  // Setup the stepper ----------------------------------------
112  auto model = rcp(new Tempus_Test::SinCosModel<double>());
113  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double> >(model);
114  auto stepper = rcp(new Tempus::StepperSubcycling<double>());
115  auto stepperBE = Tempus::createStepperBackwardEuler(modelME, Teuchos::null);
116  stepper->setSubcyclingStepper(stepperBE);
117 
118  // Setup SolutionHistory ------------------------------------
119  auto inArgsIC = model->getNominalValues();
120  auto icSolution =
121  rcp_const_cast<Thyra::VectorBase<double> >(inArgsIC.get_x());
122  auto icState = Tempus::createSolutionStateX(icSolution);
123  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
124  solutionHistory->addState(icState);
125  solutionHistory->initWorkingState();
126 
127  // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
128  stepper->setInitialConditions(solutionHistory);
129  stepper->initialize();
130 
131  // Test
132  stepper->setSubcyclingInitTimeStep(0.25);
133  stepper->setSubcyclingMaxTimeStep(0.5);
134  double maxTimeStep_Set = stepper->getSubcyclingMaxTimeStep();
135  stepper->takeStep(solutionHistory);
136  double maxTimeStep_After = stepper->getSubcyclingMaxTimeStep();
137 
138  TEST_FLOATING_EQUALITY(maxTimeStep_Set, maxTimeStep_After, 1.0e-14);
139 }
140 
141 // ************************************************************
142 // ************************************************************
143 class StepperSubcyclingModifierTest
144  : virtual public Tempus::StepperSubcyclingModifierBase<double> {
145  public:
147  StepperSubcyclingModifierTest()
148  : testBEGIN_STEP(false),
149  testEND_STEP(false),
150  testCurrentValue(-0.99),
151  testWorkingValue(-0.99),
152  testDt(1.5),
153  testName("")
154  {
155  }
157  virtual ~StepperSubcyclingModifierTest() {}
158 
160  virtual void modify(
164  actLoc)
165  {
166  switch (actLoc) {
168  testBEGIN_STEP = true;
169  auto x = sh->getCurrentState()->getX();
170  testCurrentValue = get_ele(*(x), 0);
171  testName = "Subcycling - Modifier";
172  stepper->setStepperName(testName);
173  break;
174  }
176  testEND_STEP = true;
177  auto x = sh->getWorkingState()->getX();
178  testWorkingValue = get_ele(*(x), 0);
179  testDt = sh->getWorkingState()->getTimeStep() / 10.0;
180  sh->getWorkingState()->setTimeStep(testDt);
181  break;
182  }
183  default:
184  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
185  "Error - unknown action location.\n");
186  }
187  }
188 
189  bool testBEGIN_STEP;
190  bool testEND_STEP;
191  double testCurrentValue;
192  double testWorkingValue;
193  double testDt;
194  std::string testName;
195 };
196 
197 TEUCHOS_UNIT_TEST(Subcycling, AppAction_Modifier)
198 {
199  // Setup the SinCosModel ------------------------------------
200  auto model = rcp(new Tempus_Test::SinCosModel<double>());
201  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double> >(model);
202 
203  // Setup Stepper for field solve ----------------------------
204  auto stepper = rcp(new Tempus::StepperSubcycling<double>());
205  auto stepperFE = Tempus::createStepperForwardEuler(modelME, Teuchos::null);
206  auto modifier = rcp(new StepperSubcyclingModifierTest());
207  stepper->setAppAction(modifier);
208  stepper->setSubcyclingStepper(stepperFE);
209 
210  stepper->setSubcyclingMinTimeStep(15);
211  stepper->setSubcyclingInitTimeStep(15.0);
212  stepper->setSubcyclingMaxTimeStep(15.0);
213  stepper->setSubcyclingMaxFailures(10);
214  stepper->setSubcyclingMaxConsecFailures(5);
215  stepper->setSubcyclingScreenOutputIndexInterval(1);
216  stepper->setSubcyclingPrintDtChanges(true);
217 
218  // Setup TimeStepControl ------------------------------------
219  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
220  timeStepControl->setInitIndex(0);
221  timeStepControl->setInitTime(0.0);
222  timeStepControl->setFinalTime(1.0);
223  timeStepControl->setInitTimeStep(15.0);
224  timeStepControl->initialize();
225 
226  // Setup initial condition SolutionState --------------------
227  auto inArgsIC = model->getNominalValues();
228  auto icSolution =
229  rcp_const_cast<Thyra::VectorBase<double> >(inArgsIC.get_x());
230  auto icState = Tempus::createSolutionStateX(icSolution);
231  icState->setTime(timeStepControl->getInitTime());
232  icState->setIndex(timeStepControl->getInitIndex());
233  icState->setTimeStep(0.0); // dt for ICs are indicated by zero.
234  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
235 
236  // Setup SolutionHistory ------------------------------------
237  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
238  solutionHistory->setName("Forward States");
239  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
240  solutionHistory->setStorageLimit(2);
241  solutionHistory->addState(icState);
242 
243  // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
244  stepper->setInitialConditions(solutionHistory);
245  stepper->initialize();
246 
247  // Take one time step.
248  stepper->setInitialConditions(solutionHistory);
249  solutionHistory->initWorkingState();
250  solutionHistory->getWorkingState()->setTimeStep(15.0);
251  stepper->takeStep(solutionHistory);
252 
253  // Testing that each ACTION_LOCATION has been called.
254  TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
255  TEST_COMPARE(modifier->testEND_STEP, ==, true);
256 
257  // Testing that values can be set through the Modifier.
258  auto x = solutionHistory->getCurrentState()->getX();
259  TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
260  x = solutionHistory->getWorkingState()->getX();
261  TEST_FLOATING_EQUALITY(modifier->testWorkingValue, get_ele(*(x), 0), 1.0e-14);
262  auto Dt = solutionHistory->getWorkingState()->getTimeStep();
263  TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
264 
265  TEST_COMPARE(modifier->testName, ==, "Subcycling - Modifier");
266 }
267 
268 // ************************************************************
269 // ************************************************************
270 class StepperSubcyclingObserverTest
271  : virtual public Tempus::StepperSubcyclingObserverBase<double> {
272  public:
274  StepperSubcyclingObserverTest()
275  : testBEGIN_STEP(false),
276  testEND_STEP(false),
277  testCurrentValue(-0.99),
278  testWorkingValue(-0.99),
279  testDt(15.0),
280  testName("Subcyling")
281  {
282  }
283 
285  virtual ~StepperSubcyclingObserverTest() {}
286 
288  virtual void observe(
292  actLoc)
293  {
294  switch (actLoc) {
296  testBEGIN_STEP = true;
297  auto x = sh->getCurrentState()->getX();
298  testCurrentValue = get_ele(*(x), 0);
299  break;
300  }
302  testEND_STEP = true;
303  auto x = sh->getWorkingState()->getX();
304  testWorkingValue = get_ele(*(x), 0);
305  break;
306  }
307  default:
308  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
309  "Error - unknown action location.\n");
310  }
311  }
312 
313  bool testBEGIN_STEP;
314  bool testEND_STEP;
315  double testCurrentValue;
316  double testWorkingValue;
317  double testDt;
318  std::string testName;
319 };
320 
321 TEUCHOS_UNIT_TEST(Subcycling, AppAction_Observer)
322 {
323  // Setup the SinCosModel ------------------------------------
324  auto model = rcp(new Tempus_Test::SinCosModel<double>());
325  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double> >(model);
326 
327  // Setup Stepper for field solve ----------------------------
328  auto stepper = rcp(new Tempus::StepperSubcycling<double>());
329  auto stepperFE = Tempus::createStepperForwardEuler(modelME, Teuchos::null);
330  auto observer = rcp(new StepperSubcyclingObserverTest());
331  stepper->setAppAction(observer);
332  stepper->setSubcyclingStepper(stepperFE);
333 
334  stepper->setSubcyclingMinTimeStep(15);
335  stepper->setSubcyclingInitTimeStep(15.0);
336  stepper->setSubcyclingMaxTimeStep(15.0);
337  stepper->setSubcyclingMaxFailures(10);
338  stepper->setSubcyclingMaxConsecFailures(5);
339  stepper->setSubcyclingScreenOutputIndexInterval(1);
340  stepper->setSubcyclingPrintDtChanges(true);
341 
342  // Setup TimeStepControl ------------------------------------
343  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
344  timeStepControl->setInitIndex(0);
345  timeStepControl->setInitTime(0.0);
346  timeStepControl->setFinalTime(1.0);
347  timeStepControl->setInitTimeStep(15.0);
348  timeStepControl->initialize();
349 
350  // Setup initial condition SolutionState --------------------
351  auto inArgsIC = model->getNominalValues();
352  auto icSolution =
353  rcp_const_cast<Thyra::VectorBase<double> >(inArgsIC.get_x());
354  auto icState = Tempus::createSolutionStateX(icSolution);
355  icState->setTime(timeStepControl->getInitTime());
356  icState->setIndex(timeStepControl->getInitIndex());
357  icState->setTimeStep(0.0); // dt for ICs are indicated by zero.
358  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
359 
360  // Setup SolutionHistory ------------------------------------
361  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
362  solutionHistory->setName("Forward States");
363  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
364  solutionHistory->setStorageLimit(2);
365  solutionHistory->addState(icState);
366 
367  // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
368  stepper->setInitialConditions(solutionHistory);
369  stepper->initialize();
370 
371  // Take one time step.
372  stepper->setInitialConditions(solutionHistory);
373  solutionHistory->initWorkingState();
374  solutionHistory->getWorkingState()->setTimeStep(15.0);
375  stepper->takeStep(solutionHistory);
376 
377  // Testing that each ACTION_LOCATION has been called.
378  TEST_COMPARE(observer->testBEGIN_STEP, ==, true);
379  TEST_COMPARE(observer->testEND_STEP, ==, true);
380 
381  // Testing that values can be observed through the observer.
382  auto x = solutionHistory->getCurrentState()->getX();
383  TEST_FLOATING_EQUALITY(observer->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
384  x = solutionHistory->getWorkingState()->getX();
385  TEST_FLOATING_EQUALITY(observer->testWorkingValue, get_ele(*(x), 0), 1.0e-14);
386  TEST_FLOATING_EQUALITY(observer->testDt, 15.0, 1.0e-14);
387 
388  TEST_COMPARE(observer->testName, ==, "Subcyling");
389 }
390 
391 // ************************************************************
392 // ************************************************************
393 class StepperSubcyclingModifierXTest
394  : virtual public Tempus::StepperSubcyclingModifierXBase<double> {
395  public:
397  StepperSubcyclingModifierXTest()
398  : testX_BEGIN_STEP(false),
399  testXDOT_END_STEP(false),
400  testX(-0.99),
401  testXDot(-0.99),
402  testDt(1.5),
403  testTime(1.5)
404  {
405  }
406 
408  virtual ~StepperSubcyclingModifierXTest() {}
409 
412  const double time, const double dt,
414  double>::MODIFIER_TYPE modType)
415  {
416  switch (modType) {
418  testX_BEGIN_STEP = true;
419  testX = get_ele(*(x), 0);
420  testDt = dt;
421  break;
422  }
424  testXDOT_END_STEP = true;
425  testXDot = get_ele(*(x), 0);
426  testTime = time;
427  break;
428  }
429  default:
430  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
431  "Error - unknown action location.\n");
432  }
433  }
434 
435  bool testX_BEGIN_STEP;
436  bool testXDOT_END_STEP;
437  double testX;
438  double testXDot;
439  double testDt;
440  double testTime;
441 };
442 
443 TEUCHOS_UNIT_TEST(Subcycling, AppAction_ModifierX)
444 {
445  // Setup the SinCosModel ------------------------------------
446  auto model = rcp(new Tempus_Test::SinCosModel<double>());
447  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double> >(model);
448 
449  // Setup Stepper for field solve ----------------------------
450  auto stepper = rcp(new Tempus::StepperSubcycling<double>());
451  auto stepperFE = Tempus::createStepperForwardEuler(modelME, Teuchos::null);
452  auto modifierX = rcp(new StepperSubcyclingModifierXTest());
453  stepper->setAppAction(modifierX);
454  stepper->setSubcyclingStepper(stepperFE);
455 
456  stepper->setSubcyclingMinTimeStep(15);
457  stepper->setSubcyclingInitTimeStep(15.0);
458  stepper->setSubcyclingMaxTimeStep(15.0);
459  stepper->setSubcyclingMaxFailures(10);
460  stepper->setSubcyclingMaxConsecFailures(5);
461  stepper->setSubcyclingScreenOutputIndexInterval(1);
462  stepper->setSubcyclingPrintDtChanges(true);
463 
464  // Setup TimeStepControl ------------------------------------
465  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
466  timeStepControl->setInitIndex(0);
467  timeStepControl->setInitTime(0.0);
468  timeStepControl->setFinalTime(1.0);
469  timeStepControl->setInitTimeStep(15.0);
470  timeStepControl->initialize();
471 
472  // Setup initial condition SolutionState --------------------
473  auto inArgsIC = model->getNominalValues();
474  auto icSolution =
475  rcp_const_cast<Thyra::VectorBase<double> >(inArgsIC.get_x());
476  auto icSolutionDot =
477  rcp_const_cast<Thyra::VectorBase<double> >(inArgsIC.get_x_dot());
478  auto icState = Tempus::createSolutionStateX(icSolution, icSolutionDot);
479  icState->setTime(timeStepControl->getInitTime());
480  icState->setIndex(timeStepControl->getInitIndex());
481  icState->setTimeStep(0.0); // dt for ICs are indicated by zero.
482  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
483 
484  // Setup SolutionHistory ------------------------------------
485  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
486  solutionHistory->setName("Forward States");
487  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
488  solutionHistory->setStorageLimit(2);
489  solutionHistory->addState(icState);
490 
491  // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
492  stepper->setInitialConditions(solutionHistory);
493  stepper->initialize();
494 
495  // Take one time step.
496  stepper->setInitialConditions(solutionHistory);
497  solutionHistory->initWorkingState();
498  solutionHistory->getWorkingState()->setTimeStep(15.0);
499  stepper->takeStep(solutionHistory);
500 
501  // Take one time step.
502  stepper->setInitialConditions(solutionHistory);
503  solutionHistory->initWorkingState();
504  solutionHistory->getWorkingState()->setTimeStep(15.0);
505  stepper->takeStep(solutionHistory);
506 
507  // Testing that each ACTION_LOCATION has been called.
508  TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
509  TEST_COMPARE(modifierX->testXDOT_END_STEP, ==, true);
510 
511  // Testing that values can be set through the Modifier.
512  auto x = solutionHistory->getCurrentState()->getX();
513  TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
514  // Temporary memory for xDot is not guarranteed to exist outside the Stepper.
515  auto xDot = solutionHistory->getWorkingState()->getXDot();
516  if (xDot == Teuchos::null) xDot = stepper->getStepperXDot();
517 
518  TEST_FLOATING_EQUALITY(modifierX->testXDot, get_ele(*(xDot), 0), 1.0e-14);
519  auto Dt = solutionHistory->getWorkingState()->getTimeStep();
520  TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
521 
522  auto time = solutionHistory->getWorkingState()->getTime();
523  TEST_FLOATING_EQUALITY(modifierX->testTime, time, 1.0e-14);
524 }
525 
526 } // namespace Tempus_Unit_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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.
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)
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
virtual void modify(Teuchos::RCP< Thyra::VectorBase< double > >, const double, const double, const MODIFIER_TYPE modType)=0
Modify solution based on the MODIFIER_TYPE.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
ACTION_LOCATION
Indicates the location of application action (see algorithm).
Keep a fix number of states.
Teuchos::RCP< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
virtual ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const =0
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)