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