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: 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 
13 
14 #include "Tempus_TimeStepControl.hpp"
15 
16 #include "Tempus_StepperNewmarkImplicitAForm.hpp"
21 
22 #include "../TestModels/HarmonicOscillatorModel.hpp"
23 
24 namespace Tempus_Unit_Test {
25 
27 using Teuchos::RCP;
28 using Teuchos::rcp;
29 using Teuchos::rcp_const_cast;
30 using Teuchos::rcp_dynamic_cast;
31 using Teuchos::sublist;
32 
34 
35 // ************************************************************
36 // ************************************************************
37 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, Default_Construction)
38 {
40 
41  // Default construction.
43  stepper->setModel(model);
44  stepper->initialize();
45  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
46 
47  auto modifier =
49  auto modifierX =
51 
52  // Default values for construction.
53  auto solver = rcp(new Thyra::NOXNonlinearSolver());
54  solver->setParameterList(Tempus::defaultSolverParameters());
55 
56  bool useFSAL = stepper->getUseFSAL();
57  std::string ICConsistency = stepper->getICConsistency();
58  bool ICConsistencyCheck = stepper->getICConsistencyCheck();
59  bool zeroInitialGuess = stepper->getZeroInitialGuess();
60  std::string schemeName = "Average Acceleration";
61  double beta = 0.25;
62  double gamma = 0.5;
63 
64  // Test the set functions.
65  stepper->setAppAction(modifier);
66  stepper->initialize();
67  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
68  stepper->setAppAction(modifierX);
69  stepper->initialize();
70  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
71  stepper->setSolver(solver);
72  stepper->initialize();
73  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
74  stepper->setUseFSAL(useFSAL);
75  stepper->initialize();
76  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
77  stepper->setICConsistency(ICConsistency);
78  stepper->initialize();
79  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
80  stepper->setICConsistencyCheck(ICConsistencyCheck);
81  stepper->initialize();
82  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
83  stepper->setZeroInitialGuess(zeroInitialGuess);
84  stepper->initialize();
85  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
86 
87  stepper->setSchemeName(schemeName);
88  stepper->initialize();
89  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
90  stepper->setBeta(beta);
91  stepper->initialize();
92  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
93  stepper->setGamma(gamma);
94  stepper->initialize();
95  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
96 
97  // Full argument list construction.
99  model, solver, useFSAL, ICConsistency, ICConsistencyCheck,
100  zeroInitialGuess, schemeName, beta, gamma, modifier));
101  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
102 
103  // Test stepper properties.
104  TEUCHOS_ASSERT(stepper->getOrder() == 2);
105 }
106 
107 // ************************************************************
108 // ************************************************************
109 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, StepperFactory_Construction)
110 {
112  testFactoryConstruction("Newmark Implicit a-Form", model);
113 }
114 
115 // ************************************************************
116 // ************************************************************
117 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, ConstructingFromDefaults)
118 {
119  double dt = 1.0;
120  std::vector<std::string> options;
121  options.push_back("Default Parameters");
122  options.push_back("ICConsistency and Check");
123 
124  for (const auto& option : options) {
125  // Read params from .xml file
126  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
127  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder."
128  "xml");
129  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
130 
131  // Setup the HarmonicOscillatorModel
132  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
133  auto model =
135  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
136 
137  // Setup Stepper for field solve ----------------------------
139  Tempus::createStepperNewmarkImplicitAForm(modelME, Teuchos::null);
140  if (option == "ICConsistency and Check") {
141  stepper->setICConsistency("Consistent");
142  stepper->setICConsistencyCheck(true);
143  }
144  stepper->initialize();
145 
146  // Setup TimeStepControl ------------------------------------
147  RCP<Tempus::TimeStepControl<double>> timeStepControl =
149  ParameterList tscPL =
150  pl->sublist("Default Integrator").sublist("Time Step Control");
151  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
152  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
153  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
154  timeStepControl->setInitTimeStep(dt);
155  timeStepControl->initialize();
156 
157  // Setup initial condition SolutionState --------------------
158  using Teuchos::rcp_const_cast;
159  auto inArgsIC = model->getNominalValues();
161  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
163  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
164  RCP<Thyra::VectorBase<double>> icXDotDot =
165  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
167  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
168  icState->setTime(timeStepControl->getInitTime());
169  icState->setIndex(timeStepControl->getInitIndex());
170  icState->setTimeStep(0.0);
171  icState->setOrder(stepper->getOrder());
172  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
173 
174  // Setup SolutionHistory ------------------------------------
175  RCP<Tempus::SolutionHistory<double>> solutionHistory =
177  solutionHistory->setName("Forward States");
178  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
179  solutionHistory->setStorageLimit(2);
180  solutionHistory->addState(icState);
181 
182  // Ensure ICs are consistent.
183  stepper->setInitialConditions(solutionHistory);
184 
185  // Setup Integrator -----------------------------------------
187  Tempus::createIntegratorBasic<double>();
188  integrator->setStepper(stepper);
189  integrator->setTimeStepControl(timeStepControl);
190  integrator->setSolutionHistory(solutionHistory);
191  // integrator->setObserver(...);
192  integrator->initialize();
193 
194  // Integrate to timeMax
195  bool integratorStatus = integrator->advanceTime();
196  TEST_ASSERT(integratorStatus)
197 
198  // Test if at 'Final Time'
199  double time = integrator->getTime();
200  double timeFinal = pl->sublist("Default Integrator")
201  .sublist("Time Step Control")
202  .get<double>("Final Time");
203  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
204 
205  // Time-integrated solution and the exact solution
206  RCP<Thyra::VectorBase<double>> x = integrator->getX();
208  model->getExactSolution(time).get_x();
209 
210  // Calculate the error
211  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
212  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
213 
214  // Check the order and intercept
215  out << " Stepper = " << stepper->description() << "\n with "
216  << option << std::endl;
217  out << " =========================" << std::endl;
218  out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
219  out << " Computed solution: " << get_ele(*(x), 0) << std::endl;
220  out << " Difference : " << get_ele(*(xdiff), 0) << std::endl;
221  out << " =========================" << std::endl;
222  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4);
223  }
224 }
225 
226 // ************************************************************
227 // ************************************************************
228 class StepperNewmarkImplicitAFormModifierTest
229  : virtual public Tempus::StepperNewmarkImplicitAFormModifierBase<double> {
230  public:
232  StepperNewmarkImplicitAFormModifierTest()
233  : testBEGIN_STEP(false),
234  testBEFORE_SOLVE(false),
235  testAFTER_SOLVE(false),
236  testEND_STEP(false),
237  testCurrentValue(-0.99),
238  testDt(-1.5),
239  testName("")
240  {
241  }
242 
244  virtual ~StepperNewmarkImplicitAFormModifierTest() {}
245 
247  virtual void modify(
251  double>::ACTION_LOCATION actLoc)
252  {
253  switch (actLoc) {
255  testBEGIN_STEP = true;
256  break;
257  }
259  testBEFORE_SOLVE = true;
260  testDt = sh->getWorkingState()->getTimeStep();
261  // sh->getWorkingState()->setTimeStep(testDt);
262  break;
263  }
265  testAFTER_SOLVE = true;
266  testName = "Newmark Implicit A Form - Modifier";
267  stepper->setStepperName(testName);
268  break;
269  }
271  testEND_STEP = true;
272  auto x = sh->getWorkingState()->getX();
273  testCurrentValue = get_ele(*(x), 0);
274  break;
275  }
276  default:
277  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
278  "Error - unknown action location.\n");
279  }
280  }
281 
282  bool testBEGIN_STEP;
283  bool testBEFORE_SOLVE;
284  bool testAFTER_SOLVE;
285  bool testEND_STEP;
286  double testCurrentValue;
287  double testDt;
288  std::string testName;
289 };
290 
291 // ************************************************************
292 // ************************************************************
293 class StepperNewmarkImplicitAFormModifierXTest
294  : virtual public Tempus::StepperNewmarkImplicitAFormModifierXBase<double> {
295  public:
297  StepperNewmarkImplicitAFormModifierXTest()
298  : testX_BEGIN_STEP(false),
299  testX_BEFORE_SOLVE(false),
300  testX_AFTER_SOLVE(false),
301  testX_END_STEP(false),
302  testX(-0.99),
303  testXDot(-0.99),
304  testDt(-1.5),
305  testTime(-1.5)
306  {
307  }
308 
310  virtual ~StepperNewmarkImplicitAFormModifierXTest() {}
311 
313  virtual void modify(
314  Teuchos::RCP<Thyra::VectorBase<double>> x, const double time,
315  const double dt,
317  double>::MODIFIER_TYPE modType)
318  {
319  switch (modType) {
321  testX_BEGIN_STEP = true;
322  testX = get_ele(*(x), 0);
323  break;
324  }
326  testX_BEFORE_SOLVE = true;
327  testDt = dt;
328  break;
329  }
331  testX_AFTER_SOLVE = true;
332  testTime = time;
333  testX = get_ele(*(x), 0);
334  break;
335  }
337  testX_END_STEP = true;
338  testX = get_ele(*(x), 0);
339  break;
340  }
341  default:
342  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
343  "Error - unknown action location.\n");
344  }
345  }
346 
347  bool testX_BEGIN_STEP;
348  bool testX_BEFORE_SOLVE;
349  bool testX_AFTER_SOLVE;
350  bool testX_END_STEP;
351  double testX;
352  double testXDot;
353  double testDt;
354  double testTime;
355 };
356 
357 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, AppAction_Modifier)
358 {
360  using Teuchos::RCP;
361  using Teuchos::sublist;
362 
363  double dt = 1.0;
364 
365  // Read params from .xml file
366  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
367  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
368  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
369 
370  // Setup the HarmonicOscillatorModel
371  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
374 
375  // Setup Stepper for field solve ----------------------------
377  Tempus::createStepperNewmarkImplicitAForm(model, Teuchos::null);
378 
379  auto modifier = rcp(new StepperNewmarkImplicitAFormModifierTest());
380  stepper->setAppAction(modifier);
381  stepper->initialize();
382 
383  // Setup TimeStepControl ------------------------------------
384  RCP<Tempus::TimeStepControl<double>> timeStepControl =
386  ParameterList tscPL =
387  pl->sublist("Default Integrator").sublist("Time Step Control");
388  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
389  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
390  timeStepControl->setFinalTime(dt);
391  timeStepControl->setInitTimeStep(dt);
392  timeStepControl->initialize();
393 
394  // Setup initial condition SolutionState --------------------
395  using Teuchos::rcp_const_cast;
396  auto inArgsIC = model->getNominalValues();
398  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
400  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
401  RCP<Thyra::VectorBase<double>> icXDotDot =
402  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
404  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
405  icState->setTime(timeStepControl->getInitTime());
406  icState->setIndex(timeStepControl->getInitIndex());
407  icState->setTimeStep(0.0);
408  icState->setOrder(stepper->getOrder());
409  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
410 
411  // Setup SolutionHistory ------------------------------------
412  RCP<Tempus::SolutionHistory<double>> solutionHistory =
414  solutionHistory->setName("Forward States");
415  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
416  solutionHistory->setStorageLimit(2);
417  solutionHistory->addState(icState);
418 
419  // Setup Integrator -----------------------------------------
421  Tempus::createIntegratorBasic<double>();
422  integrator->setStepper(stepper);
423  integrator->setTimeStepControl(timeStepControl);
424  integrator->setSolutionHistory(solutionHistory);
425  integrator->initialize();
426 
427  // Integrate to timeMax
428  bool integratorStatus = integrator->advanceTime();
429  TEST_ASSERT(integratorStatus)
430 
431  // Testing that each ACTION_LOCATION has been called.
432  TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
433  TEST_COMPARE(modifier->testBEFORE_SOLVE, ==, true);
434  TEST_COMPARE(modifier->testAFTER_SOLVE, ==, true);
435  TEST_COMPARE(modifier->testEND_STEP, ==, true);
436 
437  // Testing that values can be set through the Modifier.
438  auto x = integrator->getX();
439  auto Dt = integrator->getTime();
440  TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
441  TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
442  TEST_COMPARE(modifier->testName, ==, stepper->getStepperName());
443 }
444 
445 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, AppAction_ModifierX)
446 {
448  using Teuchos::RCP;
449  using Teuchos::sublist;
450 
451  double dt = 1.0;
452 
453  // Read params from .xml file
454  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
455  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
456  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
457 
458  // Setup the HarmonicOscillatorModel
459  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
462 
463  // Setup Stepper for field solve ----------------------------
465  Tempus::createStepperNewmarkImplicitAForm(model, Teuchos::null);
466 
467  auto modifierX = rcp(new StepperNewmarkImplicitAFormModifierXTest());
468  stepper->setAppAction(modifierX);
469  stepper->initialize();
470 
471  // Setup TimeStepControl ------------------------------------
472  RCP<Tempus::TimeStepControl<double>> timeStepControl =
474  ParameterList tscPL =
475  pl->sublist("Default Integrator").sublist("Time Step Control");
476  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
477  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
478  timeStepControl->setFinalTime(dt);
479  timeStepControl->setInitTimeStep(dt);
480  timeStepControl->initialize();
481 
482  // Setup initial condition SolutionState --------------------
483  using Teuchos::rcp_const_cast;
484  auto inArgsIC = model->getNominalValues();
486  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
488  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
489  RCP<Thyra::VectorBase<double>> icXDotDot =
490  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
492  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
493  icState->setTime(timeStepControl->getInitTime());
494  icState->setIndex(timeStepControl->getInitIndex());
495  icState->setTimeStep(0.0);
496  icState->setOrder(stepper->getOrder());
497  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
498 
499  // Setup SolutionHistory ------------------------------------
500  RCP<Tempus::SolutionHistory<double>> solutionHistory =
502  solutionHistory->setName("Forward States");
503  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
504  solutionHistory->setStorageLimit(2);
505  solutionHistory->addState(icState);
506 
507  // Setup Integrator -----------------------------------------
509  Tempus::createIntegratorBasic<double>();
510  integrator->setStepper(stepper);
511  integrator->setTimeStepControl(timeStepControl);
512  integrator->setSolutionHistory(solutionHistory);
513  integrator->initialize();
514 
515  // Integrate to timeMax
516  bool integratorStatus = integrator->advanceTime();
517  TEST_ASSERT(integratorStatus)
518 
519  // Testing that each ACTION_LOCATION has been called.
520  TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
521  TEST_COMPARE(modifierX->testX_BEFORE_SOLVE, ==, true);
522  TEST_COMPARE(modifierX->testX_AFTER_SOLVE, ==, true);
523  TEST_COMPARE(modifierX->testX_END_STEP, ==, true);
524 
525  // Testing that values can be set through the Modifier.
526  auto Dt = integrator->getTime();
527  TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
528 
529  const auto x = integrator->getX();
530  TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
531 }
532 
533 } // 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)