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