Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_UnitTest_NewmarkImplicitAForm.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 
13 #include "Tempus_TimeStepControl.hpp"
14 
15 #include "Tempus_StepperNewmarkImplicitAForm.hpp"
20 
21 #include "../TestModels/HarmonicOscillatorModel.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(NewmarkImplicitAForm, Default_Construction)
37 {
39 
40  // Default construction.
42  stepper->setModel(model);
43  stepper->initialize();
44  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
45 
46  auto modifier =
48  auto modifierX =
50 
51  // Default values for construction.
52  auto solver = rcp(new Thyra::NOXNonlinearSolver());
53  solver->setParameterList(Tempus::defaultSolverParameters());
54 
55  bool useFSAL = stepper->getUseFSAL();
56  std::string ICConsistency = stepper->getICConsistency();
57  bool ICConsistencyCheck = stepper->getICConsistencyCheck();
58  bool zeroInitialGuess = stepper->getZeroInitialGuess();
59  std::string schemeName = "Average Acceleration";
60  double beta = 0.25;
61  double gamma = 0.5;
62 
63  // Test the set functions.
64  stepper->setAppAction(modifier);
65  stepper->initialize();
66  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
67  stepper->setAppAction(modifierX);
68  stepper->initialize();
69  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
70  stepper->setSolver(solver);
71  stepper->initialize();
72  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
73  stepper->setUseFSAL(useFSAL);
74  stepper->initialize();
75  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
76  stepper->setICConsistency(ICConsistency);
77  stepper->initialize();
78  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
79  stepper->setICConsistencyCheck(ICConsistencyCheck);
80  stepper->initialize();
81  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
82  stepper->setZeroInitialGuess(zeroInitialGuess);
83  stepper->initialize();
84  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
85 
86  stepper->setSchemeName(schemeName);
87  stepper->initialize();
88  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
89  stepper->setBeta(beta);
90  stepper->initialize();
91  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
92  stepper->setGamma(gamma);
93  stepper->initialize();
94  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
95 
96  // Full argument list construction.
98  model, solver, useFSAL, ICConsistency, ICConsistencyCheck,
99  zeroInitialGuess, schemeName, beta, gamma, modifier));
100  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
101 
102  // Test stepper properties.
103  TEUCHOS_ASSERT(stepper->getOrder() == 2);
104 }
105 
106 // ************************************************************
107 // ************************************************************
108 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, StepperFactory_Construction)
109 {
111  testFactoryConstruction("Newmark Implicit a-Form", model);
112 }
113 
114 // ************************************************************
115 // ************************************************************
116 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, ConstructingFromDefaults)
117 {
118  double dt = 1.0;
119  std::vector<std::string> options;
120  options.push_back("Default Parameters");
121  options.push_back("ICConsistency and Check");
122 
123  for (const auto& option : options) {
124  // Read params from .xml file
125  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
126  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder."
127  "xml");
128  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
129 
130  // Setup the HarmonicOscillatorModel
131  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
132  auto model =
134  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
135 
136  // Setup Stepper for field solve ----------------------------
138  Tempus::createStepperNewmarkImplicitAForm(modelME, Teuchos::null);
139  if (option == "ICConsistency and Check") {
140  stepper->setICConsistency("Consistent");
141  stepper->setICConsistencyCheck(true);
142  }
143  stepper->initialize();
144 
145  // Setup TimeStepControl ------------------------------------
146  RCP<Tempus::TimeStepControl<double>> timeStepControl =
148  ParameterList tscPL =
149  pl->sublist("Default Integrator").sublist("Time Step Control");
150  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
151  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
152  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
153  timeStepControl->setInitTimeStep(dt);
154  timeStepControl->initialize();
155 
156  // Setup initial condition SolutionState --------------------
157  using Teuchos::rcp_const_cast;
158  auto inArgsIC = model->getNominalValues();
160  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
162  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
163  RCP<Thyra::VectorBase<double>> icXDotDot =
164  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
166  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
167  icState->setTime(timeStepControl->getInitTime());
168  icState->setIndex(timeStepControl->getInitIndex());
169  icState->setTimeStep(0.0);
170  icState->setOrder(stepper->getOrder());
171  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
172 
173  // Setup SolutionHistory ------------------------------------
174  RCP<Tempus::SolutionHistory<double>> solutionHistory =
176  solutionHistory->setName("Forward States");
177  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
178  solutionHistory->setStorageLimit(2);
179  solutionHistory->addState(icState);
180 
181  // Ensure ICs are consistent.
182  stepper->setInitialConditions(solutionHistory);
183 
184  // Setup Integrator -----------------------------------------
186  Tempus::createIntegratorBasic<double>();
187  integrator->setStepper(stepper);
188  integrator->setTimeStepControl(timeStepControl);
189  integrator->setSolutionHistory(solutionHistory);
190  // integrator->setObserver(...);
191  integrator->initialize();
192 
193  // Integrate to timeMax
194  bool integratorStatus = integrator->advanceTime();
195  TEST_ASSERT(integratorStatus)
196 
197  // Test if at 'Final Time'
198  double time = integrator->getTime();
199  double timeFinal = pl->sublist("Default Integrator")
200  .sublist("Time Step Control")
201  .get<double>("Final Time");
202  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
203 
204  // Time-integrated solution and the exact solution
205  RCP<Thyra::VectorBase<double>> x = integrator->getX();
207  model->getExactSolution(time).get_x();
208 
209  // Calculate the error
210  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
211  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
212 
213  // Check the order and intercept
214  out << " Stepper = " << stepper->description() << "\n with "
215  << option << std::endl;
216  out << " =========================" << std::endl;
217  out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
218  out << " Computed solution: " << get_ele(*(x), 0) << std::endl;
219  out << " Difference : " << get_ele(*(xdiff), 0) << std::endl;
220  out << " =========================" << std::endl;
221  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4);
222  }
223 }
224 
225 // ************************************************************
226 // ************************************************************
227 class StepperNewmarkImplicitAFormModifierTest
228  : virtual public Tempus::StepperNewmarkImplicitAFormModifierBase<double> {
229  public:
231  StepperNewmarkImplicitAFormModifierTest()
232  : testBEGIN_STEP(false),
233  testBEFORE_SOLVE(false),
234  testAFTER_SOLVE(false),
235  testEND_STEP(false),
236  testCurrentValue(-0.99),
237  testDt(-1.5),
238  testName("")
239  {
240  }
241 
243  virtual ~StepperNewmarkImplicitAFormModifierTest() {}
244 
246  virtual void modify(
250  double>::ACTION_LOCATION actLoc)
251  {
252  switch (actLoc) {
254  testBEGIN_STEP = true;
255  break;
256  }
258  testBEFORE_SOLVE = true;
259  testDt = sh->getWorkingState()->getTimeStep();
260  // sh->getWorkingState()->setTimeStep(testDt);
261  break;
262  }
264  testAFTER_SOLVE = true;
265  testName = "Newmark Implicit A Form - Modifier";
266  stepper->setStepperName(testName);
267  break;
268  }
270  testEND_STEP = true;
271  auto x = sh->getWorkingState()->getX();
272  testCurrentValue = get_ele(*(x), 0);
273  break;
274  }
275  default:
276  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
277  "Error - unknown action location.\n");
278  }
279  }
280 
281  bool testBEGIN_STEP;
282  bool testBEFORE_SOLVE;
283  bool testAFTER_SOLVE;
284  bool testEND_STEP;
285  double testCurrentValue;
286  double testDt;
287  std::string testName;
288 };
289 
290 // ************************************************************
291 // ************************************************************
292 class StepperNewmarkImplicitAFormModifierXTest
293  : virtual public Tempus::StepperNewmarkImplicitAFormModifierXBase<double> {
294  public:
296  StepperNewmarkImplicitAFormModifierXTest()
297  : testX_BEGIN_STEP(false),
298  testX_BEFORE_SOLVE(false),
299  testX_AFTER_SOLVE(false),
300  testX_END_STEP(false),
301  testX(-0.99),
302  testXDot(-0.99),
303  testDt(-1.5),
304  testTime(-1.5)
305  {
306  }
307 
309  virtual ~StepperNewmarkImplicitAFormModifierXTest() {}
310 
312  virtual void modify(
313  Teuchos::RCP<Thyra::VectorBase<double>> x, const double time,
314  const double dt,
316  double>::MODIFIER_TYPE modType)
317  {
318  switch (modType) {
320  testX_BEGIN_STEP = true;
321  testX = get_ele(*(x), 0);
322  break;
323  }
325  testX_BEFORE_SOLVE = true;
326  testDt = dt;
327  break;
328  }
330  testX_AFTER_SOLVE = true;
331  testTime = time;
332  testX = get_ele(*(x), 0);
333  break;
334  }
336  testX_END_STEP = true;
337  testX = get_ele(*(x), 0);
338  break;
339  }
340  default:
341  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
342  "Error - unknown action location.\n");
343  }
344  }
345 
346  bool testX_BEGIN_STEP;
347  bool testX_BEFORE_SOLVE;
348  bool testX_AFTER_SOLVE;
349  bool testX_END_STEP;
350  double testX;
351  double testXDot;
352  double testDt;
353  double testTime;
354 };
355 
356 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, AppAction_Modifier)
357 {
359  using Teuchos::RCP;
360  using Teuchos::sublist;
361 
362  double dt = 1.0;
363 
364  // Read params from .xml file
365  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
366  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
367  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
368 
369  // Setup the HarmonicOscillatorModel
370  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
373 
374  // Setup Stepper for field solve ----------------------------
376  Tempus::createStepperNewmarkImplicitAForm(model, Teuchos::null);
377 
378  auto modifier = rcp(new StepperNewmarkImplicitAFormModifierTest());
379  stepper->setAppAction(modifier);
380  stepper->initialize();
381 
382  // Setup TimeStepControl ------------------------------------
383  RCP<Tempus::TimeStepControl<double>> timeStepControl =
385  ParameterList tscPL =
386  pl->sublist("Default Integrator").sublist("Time Step Control");
387  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
388  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
389  timeStepControl->setFinalTime(dt);
390  timeStepControl->setInitTimeStep(dt);
391  timeStepControl->initialize();
392 
393  // Setup initial condition SolutionState --------------------
394  using Teuchos::rcp_const_cast;
395  auto inArgsIC = model->getNominalValues();
397  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
399  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
400  RCP<Thyra::VectorBase<double>> icXDotDot =
401  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
403  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
404  icState->setTime(timeStepControl->getInitTime());
405  icState->setIndex(timeStepControl->getInitIndex());
406  icState->setTimeStep(0.0);
407  icState->setOrder(stepper->getOrder());
408  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
409 
410  // Setup SolutionHistory ------------------------------------
411  RCP<Tempus::SolutionHistory<double>> solutionHistory =
413  solutionHistory->setName("Forward States");
414  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
415  solutionHistory->setStorageLimit(2);
416  solutionHistory->addState(icState);
417 
418  // Setup Integrator -----------------------------------------
420  Tempus::createIntegratorBasic<double>();
421  integrator->setStepper(stepper);
422  integrator->setTimeStepControl(timeStepControl);
423  integrator->setSolutionHistory(solutionHistory);
424  integrator->initialize();
425 
426  // Integrate to timeMax
427  bool integratorStatus = integrator->advanceTime();
428  TEST_ASSERT(integratorStatus)
429 
430  // Testing that each ACTION_LOCATION has been called.
431  TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
432  TEST_COMPARE(modifier->testBEFORE_SOLVE, ==, true);
433  TEST_COMPARE(modifier->testAFTER_SOLVE, ==, true);
434  TEST_COMPARE(modifier->testEND_STEP, ==, true);
435 
436  // Testing that values can be set through the Modifier.
437  auto x = integrator->getX();
438  auto Dt = integrator->getTime();
439  TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
440  TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
441  TEST_COMPARE(modifier->testName, ==, stepper->getStepperName());
442 }
443 
444 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, AppAction_ModifierX)
445 {
447  using Teuchos::RCP;
448  using Teuchos::sublist;
449 
450  double dt = 1.0;
451 
452  // Read params from .xml file
453  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
454  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
455  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
456 
457  // Setup the HarmonicOscillatorModel
458  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
461 
462  // Setup Stepper for field solve ----------------------------
464  Tempus::createStepperNewmarkImplicitAForm(model, Teuchos::null);
465 
466  auto modifierX = rcp(new StepperNewmarkImplicitAFormModifierXTest());
467  stepper->setAppAction(modifierX);
468  stepper->initialize();
469 
470  // Setup TimeStepControl ------------------------------------
471  RCP<Tempus::TimeStepControl<double>> timeStepControl =
473  ParameterList tscPL =
474  pl->sublist("Default Integrator").sublist("Time Step Control");
475  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
476  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
477  timeStepControl->setFinalTime(dt);
478  timeStepControl->setInitTimeStep(dt);
479  timeStepControl->initialize();
480 
481  // Setup initial condition SolutionState --------------------
482  using Teuchos::rcp_const_cast;
483  auto inArgsIC = model->getNominalValues();
485  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
487  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
488  RCP<Thyra::VectorBase<double>> icXDotDot =
489  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
491  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
492  icState->setTime(timeStepControl->getInitTime());
493  icState->setIndex(timeStepControl->getInitIndex());
494  icState->setTimeStep(0.0);
495  icState->setOrder(stepper->getOrder());
496  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
497 
498  // Setup SolutionHistory ------------------------------------
499  RCP<Tempus::SolutionHistory<double>> solutionHistory =
501  solutionHistory->setName("Forward States");
502  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
503  solutionHistory->setStorageLimit(2);
504  solutionHistory->addState(icState);
505 
506  // Setup Integrator -----------------------------------------
508  Tempus::createIntegratorBasic<double>();
509  integrator->setStepper(stepper);
510  integrator->setTimeStepControl(timeStepControl);
511  integrator->setSolutionHistory(solutionHistory);
512  integrator->initialize();
513 
514  // Integrate to timeMax
515  bool integratorStatus = integrator->advanceTime();
516  TEST_ASSERT(integratorStatus)
517 
518  // Testing that each ACTION_LOCATION has been called.
519  TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
520  TEST_COMPARE(modifierX->testX_BEFORE_SOLVE, ==, true);
521  TEST_COMPARE(modifierX->testX_AFTER_SOLVE, ==, true);
522  TEST_COMPARE(modifierX->testX_END_STEP, ==, true);
523 
524  // Testing that values can be set through the Modifier.
525  auto Dt = integrator->getTime();
526  TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
527 
528  const auto x = integrator->getX();
529  TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
530 }
531 
532 } // 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.
Application Action for StepperNewmarkImplicitAForm.
Teuchos::RCP< StepperNewmarkImplicitAForm< Scalar > > createStepperNewmarkImplicitAForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
void testFactoryConstruction(std::string stepperType, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model)
Unit test utility for Stepper construction through StepperFactory.
T & get(const std::string &name, T def_value)
#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)
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< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.
Consider the ODE: where is a constant, is a constant damping parameter, is a constant forcing par...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ptr< T > ptr() const
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
virtual void modify(Teuchos::RCP< SolutionHistory< double > >, Teuchos::RCP< StepperNewmarkImplicitAForm< double > >, const typename StepperNewmarkImplicitAFormAppAction< double >::ACTION_LOCATION actLoc)=0
Modify NewmarkImplicitAForm Stepper.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEUCHOS_ASSERT(assertion_test)
Newmark time stepper in acceleration form (a-form).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)