Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_UnitTest_NewmarkImplicitDForm.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_StepperNewmarkImplicitDForm.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(NewmarkImplicitDForm, Constructing_From_Defaults)
38 {
39  double dt = 1.0;
40 
41  // Read params from .xml file
42  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
43  "Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.xml");
44  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
45 
46  // Setup the HarmonicOscillatorModel
47  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
48  auto model = Teuchos::rcp(
50  auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
51 
52  // Setup Stepper for field solve ----------------------------
53  auto stepper =
54  Tempus::createStepperNewmarkImplicitDForm(modelME, Teuchos::null);
55 
56  // Setup TimeStepControl ------------------------------------
57  RCP<Tempus::TimeStepControl<double>> timeStepControl =
59  ParameterList tscPL =
60  pl->sublist("Default Integrator").sublist("Time Step Control");
61  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
62  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
63  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
64  timeStepControl->setInitTimeStep(dt);
65  timeStepControl->initialize();
66 
67  // Setup initial condition SolutionState --------------------
68  using Teuchos::rcp_const_cast;
69  auto inArgsIC = model->getNominalValues();
71  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
73  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
75  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
77  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
78  icState->setTime(timeStepControl->getInitTime());
79  icState->setIndex(timeStepControl->getInitIndex());
80  icState->setTimeStep(0.0);
81  icState->setOrder(stepper->getOrder());
82  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
83 
84  // Setup SolutionHistory ------------------------------------
85  RCP<Tempus::SolutionHistory<double>> solutionHistory =
87  solutionHistory->setName("Forward States");
88  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
89  solutionHistory->setStorageLimit(2);
90  solutionHistory->addState(icState);
91 
92  // Setup Integrator -----------------------------------------
94  Tempus::createIntegratorBasic<double>();
95  integrator->setStepper(stepper);
96  integrator->setTimeStepControl(timeStepControl);
97  integrator->setSolutionHistory(solutionHistory);
98  // integrator->setObserver(...);
99  integrator->initialize();
100 
101  // Integrate to timeMax
102  bool integratorStatus = integrator->advanceTime();
103  TEST_ASSERT(integratorStatus)
104 
105  // Test if at 'Final Time'
106  double time = integrator->getTime();
107  double timeFinal = pl->sublist("Default Integrator")
108  .sublist("Time Step Control")
109  .get<double>("Final Time");
110  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
111 
112  // Time-integrated solution and the exact solution
113  RCP<Thyra::VectorBase<double>> x = integrator->getX();
115  model->getExactSolution(time).get_x();
116 
117  // Calculate the error
118  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
119  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
120 
121  // Check the order and intercept
122  out << " Stepper = " << stepper->description() << std::endl;
123  out << " =========================" << std::endl;
124  out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
125  out << " Computed solution: " << get_ele(*(x), 0) << std::endl;
126  out << " Difference : " << get_ele(*(xdiff), 0) << std::endl;
127  out << " =========================" << std::endl;
128  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4);
129 }
130 
131 // ************************************************************
132 class StepperNewmarkImplicitDFormModifierTest
133  : virtual public Tempus::StepperNewmarkImplicitDFormModifierBase<double> {
134  public:
136  StepperNewmarkImplicitDFormModifierTest()
137  : testBEGIN_STEP(false),
138  testBEFORE_SOLVE(false),
139  testAFTER_SOLVE(false),
140  testEND_STEP(false),
141  testCurrentValue(-0.99),
142  testDt(-1.5),
143  testName("")
144  {
145  }
146 
148  virtual ~StepperNewmarkImplicitDFormModifierTest() {}
149 
151  virtual void modify(
155  double>::ACTION_LOCATION actLoc)
156  {
157  switch (actLoc) {
159  testBEGIN_STEP = true;
160  break;
161  }
163  testBEFORE_SOLVE = true;
164  testDt = sh->getWorkingState()->getTimeStep();
165  // sh->getWorkingState()->setTimeStep(testDt);
166  break;
167  }
169  testAFTER_SOLVE = true;
170  testName = "Newmark Implicit A Form - Modifier";
171  stepper->setStepperName(testName);
172  break;
173  }
175  testEND_STEP = true;
176  auto x = sh->getWorkingState()->getX();
177  testCurrentValue = get_ele(*(x), 0);
178  break;
179  }
180  default:
181  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
182  "Error - unknown action location.\n");
183  }
184  }
185 
186  bool testBEGIN_STEP;
187  bool testBEFORE_SOLVE;
188  bool testAFTER_SOLVE;
189  bool testEND_STEP;
190  double testCurrentValue;
191  double testDt;
192  std::string testName;
193 };
194 
195 // ************************************************************
196 // ************************************************************
197 class StepperNewmarkImplicitDFormModifierXTest
198  : virtual public Tempus::StepperNewmarkImplicitDFormModifierXBase<double> {
199  public:
201  StepperNewmarkImplicitDFormModifierXTest()
202  : testX_BEGIN_STEP(false),
203  testX_BEFORE_SOLVE(false),
204  testX_AFTER_SOLVE(false),
205  testX_END_STEP(false),
206  testX(-0.99),
207  testXDot(-0.99),
208  testDt(-1.5),
209  testTime(-1.5)
210  {
211  }
212 
214  virtual ~StepperNewmarkImplicitDFormModifierXTest() {}
215 
217  virtual void modify(
218  Teuchos::RCP<Thyra::VectorBase<double>> x, const double time,
219  const double dt,
221  double>::MODIFIER_TYPE modType)
222  {
223  switch (modType) {
225  testX_BEGIN_STEP = true;
226  testX = get_ele(*(x), 0);
227  break;
228  }
230  testX_BEFORE_SOLVE = true;
231  testDt = dt;
232  break;
233  }
235  testX_AFTER_SOLVE = true;
236  testTime = time;
237  testX = get_ele(*(x), 0);
238  break;
239  }
241  testX_END_STEP = true;
242  testX = get_ele(*(x), 0);
243  break;
244  }
245  default:
246  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
247  "Error - unknown action location.\n");
248  }
249  }
250 
251  bool testX_BEGIN_STEP;
252  bool testX_BEFORE_SOLVE;
253  bool testX_AFTER_SOLVE;
254  bool testX_END_STEP;
255  double testX;
256  double testXDot;
257  double testDt;
258  double testTime;
259 };
260 
261 // ************************************************************
262 // ************************************************************
263 TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, Default_Construction)
264 {
266 
267  // Default construction.
268  auto stepper = rcp(new Tempus::StepperNewmarkImplicitDForm<double>());
269  stepper->setModel(model);
270  stepper->initialize();
271  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
272 
273  auto modifier =
275  auto modifierX =
277 
278  // Default values for construction.
279  auto solver = rcp(new Thyra::NOXNonlinearSolver());
280  solver->setParameterList(Tempus::defaultSolverParameters());
281 
282  bool useFSAL = stepper->getUseFSAL();
283  std::string ICConsistency = stepper->getICConsistency();
284  bool ICConsistencyCheck = stepper->getICConsistencyCheck();
285  bool zeroInitialGuess = stepper->getZeroInitialGuess();
286  std::string schemeName = "Average Acceleration";
287  double beta = 0.25;
288  double gamma = 0.5;
289 
290  // Test the set functions.
291  stepper->setAppAction(modifier);
292  stepper->initialize();
293  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
294  stepper->setAppAction(modifierX);
295  stepper->initialize();
296  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
297  stepper->setSolver(solver);
298  stepper->initialize();
299  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
300  stepper->setUseFSAL(useFSAL);
301  stepper->initialize();
302  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
303  stepper->setICConsistency(ICConsistency);
304  stepper->initialize();
305  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
306  stepper->setICConsistencyCheck(ICConsistencyCheck);
307  stepper->initialize();
308  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
309  stepper->setZeroInitialGuess(zeroInitialGuess);
310  stepper->initialize();
311  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
312 
313  stepper->setSchemeName(schemeName);
314  stepper->initialize();
315  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
316  stepper->setBeta(beta);
317  stepper->initialize();
318  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
319  stepper->setGamma(gamma);
320  stepper->initialize();
321  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
322 
323  // Full argument list construction.
325  model, solver, useFSAL, ICConsistency, ICConsistencyCheck,
326  zeroInitialGuess, schemeName, beta, gamma, modifier));
327  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
328 
329  // Test stepper properties.
330  TEUCHOS_ASSERT(stepper->getOrder() == 2);
331 }
332 
333 // ************************************************************
334 // ************************************************************
335 TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, StepperFactory_Construction)
336 {
338  testFactoryConstruction("Newmark Implicit d-Form", model);
339 }
340 
341 TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, AppAction_Modifier)
342 {
344  using Teuchos::RCP;
345  using Teuchos::sublist;
346 
347  double dt = 1.0;
348 
349  // Read params from .xml file
350  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
351  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
352  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
353 
354  // Setup the HarmonicOscillatorModel
355  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
358 
359  // Setup Stepper for field solve ----------------------------
361  Tempus::createStepperNewmarkImplicitDForm(model, Teuchos::null);
362 
363  auto modifier = rcp(new StepperNewmarkImplicitDFormModifierTest());
364  stepper->setAppAction(modifier);
365  stepper->initialize();
366 
367  // Setup TimeStepControl ------------------------------------
368  RCP<Tempus::TimeStepControl<double>> timeStepControl =
370  ParameterList tscPL =
371  pl->sublist("Default Integrator").sublist("Time Step Control");
372  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
373  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
374  timeStepControl->setFinalTime(dt);
375  timeStepControl->setInitTimeStep(dt);
376  timeStepControl->initialize();
377 
378  // Setup initial condition SolutionState --------------------
379  using Teuchos::rcp_const_cast;
380  auto inArgsIC = model->getNominalValues();
382  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
384  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
385  RCP<Thyra::VectorBase<double>> icXDotDot =
386  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
388  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
389  icState->setTime(timeStepControl->getInitTime());
390  icState->setIndex(timeStepControl->getInitIndex());
391  icState->setTimeStep(0.0);
392  icState->setOrder(stepper->getOrder());
393  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
394 
395  // Setup SolutionHistory ------------------------------------
396  RCP<Tempus::SolutionHistory<double>> solutionHistory =
398  solutionHistory->setName("Forward States");
399  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
400  solutionHistory->setStorageLimit(2);
401  solutionHistory->addState(icState);
402 
403  // Setup Integrator -----------------------------------------
405  Tempus::createIntegratorBasic<double>();
406  integrator->setStepper(stepper);
407  integrator->setTimeStepControl(timeStepControl);
408  integrator->setSolutionHistory(solutionHistory);
409  integrator->initialize();
410 
411  // Integrate to timeMax
412  bool integratorStatus = integrator->advanceTime();
413  TEST_ASSERT(integratorStatus)
414 
415  // Testing that each ACTION_LOCATION has been called.
416  TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
417  TEST_COMPARE(modifier->testBEFORE_SOLVE, ==, true);
418  TEST_COMPARE(modifier->testAFTER_SOLVE, ==, true);
419  TEST_COMPARE(modifier->testEND_STEP, ==, true);
420 
421  // Testing that values can be set through the Modifier.
422  auto x = integrator->getX();
423  auto Dt = integrator->getTime();
424  TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
425  TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
426  TEST_COMPARE(modifier->testName, ==, stepper->getStepperName());
427 }
428 
429 TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, AppAction_ModifierX)
430 {
432  using Teuchos::RCP;
433  using Teuchos::sublist;
434 
435  double dt = 1.0;
436 
437  // Read params from .xml file
438  RCP<ParameterList> pList = Teuchos::getParametersFromXmlFile(
439  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
440  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
441 
442  // Setup the HarmonicOscillatorModel
443  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
446 
447  // Setup Stepper for field solve ----------------------------
449  Tempus::createStepperNewmarkImplicitDForm(model, Teuchos::null);
450 
451  auto modifierX = rcp(new StepperNewmarkImplicitDFormModifierXTest());
452  stepper->setAppAction(modifierX);
453  stepper->initialize();
454 
455  // Setup TimeStepControl ------------------------------------
456  RCP<Tempus::TimeStepControl<double>> timeStepControl =
458  ParameterList tscPL =
459  pl->sublist("Default Integrator").sublist("Time Step Control");
460  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
461  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
462  timeStepControl->setFinalTime(dt);
463  timeStepControl->setInitTimeStep(dt);
464  timeStepControl->initialize();
465 
466  // Setup initial condition SolutionState --------------------
467  using Teuchos::rcp_const_cast;
468  auto inArgsIC = model->getNominalValues();
470  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
472  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
473  RCP<Thyra::VectorBase<double>> icXDotDot =
474  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
476  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
477  icState->setTime(timeStepControl->getInitTime());
478  icState->setIndex(timeStepControl->getInitIndex());
479  icState->setTimeStep(0.0);
480  icState->setOrder(stepper->getOrder());
481  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
482 
483  // Setup SolutionHistory ------------------------------------
484  RCP<Tempus::SolutionHistory<double>> solutionHistory =
486  solutionHistory->setName("Forward States");
487  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
488  solutionHistory->setStorageLimit(2);
489  solutionHistory->addState(icState);
490 
491  // Setup Integrator -----------------------------------------
493  Tempus::createIntegratorBasic<double>();
494  integrator->setStepper(stepper);
495  integrator->setTimeStepControl(timeStepControl);
496  integrator->setSolutionHistory(solutionHistory);
497  integrator->initialize();
498 
499  // Integrate to timeMax
500  bool integratorStatus = integrator->advanceTime();
501  TEST_ASSERT(integratorStatus)
502 
503  // Testing that each ACTION_LOCATION has been called.
504  TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
505  TEST_COMPARE(modifierX->testX_BEFORE_SOLVE, ==, true);
506  TEST_COMPARE(modifierX->testX_AFTER_SOLVE, ==, true);
507  TEST_COMPARE(modifierX->testX_END_STEP, ==, true);
508 
509  // Testing that values can be set through the Modifier.
510  auto Dt = integrator->getTime();
511  TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
512 
513  const auto x = integrator->getX();
514  TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
515 }
516 
517 } // 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.
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)
Teuchos::RCP< StepperNewmarkImplicitDForm< Scalar > > createStepperNewmarkImplicitDForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar >> &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEST_COMPARE(v1, comp, v2)
virtual void modify(Teuchos::RCP< Thyra::VectorBase< double > >, const double, const double, const MODIFIER_TYPE modType)=0
Modify solution based on the MODIFIER_TYPE.
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
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...
virtual void modify(Teuchos::RCP< SolutionHistory< double > >, Teuchos::RCP< StepperNewmarkImplicitDForm< double > >, const typename StepperNewmarkImplicitDFormAppAction< double >::ACTION_LOCATION actLoc)=0
Modify NewmarkImplicitDForm Stepper.
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.
Application Action for StepperNewmarkImplicitDForm.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)