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