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