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