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 
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Thyra_VectorStdOps.hpp"
15 
16 #include "Tempus_IntegratorBasic.hpp"
17 #include "Tempus_StepperFactory.hpp"
18 #include "Tempus_SolutionHistory.hpp"
20 #include "Tempus_TimeStepControl.hpp"
21 
22 #include "Tempus_StepperNewmarkImplicitAForm.hpp"
27 
28 #include "../TestModels/SinCosModel.hpp"
29 #include "../TestModels/VanDerPolModel.hpp"
30 #include "../TestModels/HarmonicOscillatorModel.hpp"
31 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
32 
33 #include <fstream>
34 #include <vector>
35 
36 namespace Tempus_Unit_Test {
37 
38 using Teuchos::RCP;
39 using Teuchos::rcp;
40 using Teuchos::rcp_const_cast;
41 using Teuchos::rcp_dynamic_cast;
43 using Teuchos::sublist;
44 using Teuchos::getParametersFromXmlFile;
45 
47 
48 
49 // ************************************************************
50 // ************************************************************
51 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, Default_Construction)
52 {
54 
55  // Default construction.
57  stepper->setModel(model);
58  stepper->initialize();
59  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
60 
63 
64  // Default values for construction.
65  auto solver = rcp(new Thyra::NOXNonlinearSolver());
66  solver->setParameterList(Tempus::defaultSolverParameters());
67 
68  bool useFSAL = stepper->getUseFSAL();
69  std::string ICConsistency = stepper->getICConsistency();
70  bool ICConsistencyCheck = stepper->getICConsistencyCheck();
71  bool zeroInitialGuess = stepper->getZeroInitialGuess();
72  std::string schemeName = "Average Acceleration";
73  double beta = 0.25;
74  double gamma = 0.5;
75 
76 
77  // Test the set functions.
78  stepper->setAppAction(modifier); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
79  stepper->setAppAction(modifierX); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
80  stepper->setSolver(solver); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
81  stepper->setUseFSAL(useFSAL); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
82  stepper->setICConsistency(ICConsistency); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
83  stepper->setICConsistencyCheck(ICConsistencyCheck); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
84  stepper->setZeroInitialGuess(zeroInitialGuess); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
85 
86  stepper->setSchemeName(schemeName); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
87  stepper->setBeta(beta); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
88  stepper->setGamma(gamma); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
89 
90 
91  // Full argument list construction.
93  model, solver, useFSAL,
94  ICConsistency, ICConsistencyCheck, zeroInitialGuess,
95  schemeName, beta, gamma, modifier));
96  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
97 
98  // Test stepper properties.
99  TEUCHOS_ASSERT(stepper->getOrder() == 2);
100 }
101 
102 
103 // ************************************************************
104 // ************************************************************
105 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, StepperFactory_Construction)
106 {
108  testFactoryConstruction("Newmark Implicit a-Form", model);
109 }
110 
111 // ************************************************************
112 // ************************************************************
113 class StepperNewmarkImplicitAFormModifierTest
114  : virtual public Tempus::StepperNewmarkImplicitAFormModifierBase<double>
115 {
116 public:
117 
119  StepperNewmarkImplicitAFormModifierTest()
120  : testBEGIN_STEP(false), testBEFORE_SOLVE(false),
121  testAFTER_SOLVE(false), testEND_STEP(false),
122  testCurrentValue(-0.99),
123  testDt(-1.5), testName("")
124  {}
125 
127  virtual ~StepperNewmarkImplicitAFormModifierTest(){}
128 
130  virtual void modify(
134  {
135  switch(actLoc) {
137  {
138  testBEGIN_STEP = true;
139  break;
140  }
142  {
143  testBEFORE_SOLVE = true;
144  testDt = sh->getWorkingState()->getTimeStep();
145  //sh->getWorkingState()->setTimeStep(testDt);
146  break;
147  }
149  {
150  testAFTER_SOLVE = true;
151  testName = "Newmark Implicit A Form - Modifier";
152  stepper->setStepperName(testName);
153  break;
154  }
156  {
157  testEND_STEP = true;
158  auto x = sh->getWorkingState()->getX();
159  testCurrentValue = get_ele(*(x), 0);
160  break;
161  }
162  default:
163  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
164  "Error - unknown action location.\n");
165  }
166  }
167 
168  bool testBEGIN_STEP;
169  bool testBEFORE_SOLVE;
170  bool testAFTER_SOLVE;
171  bool testEND_STEP;
172  double testCurrentValue;
173  double testDt;
174  std::string testName;
175 };
176 
177 
178 // ************************************************************
179 // ************************************************************
180 class StepperNewmarkImplicitAFormModifierXTest
181  : virtual public Tempus::StepperNewmarkImplicitAFormModifierXBase<double>
182 {
183 public:
184 
186  StepperNewmarkImplicitAFormModifierXTest()
187  : testX_BEGIN_STEP(false), testX_BEFORE_SOLVE(false),
188  testX_AFTER_SOLVE(false), testX_END_STEP(false),
189  testX(-0.99), testXDot(-0.99),
190  testDt(-1.5), testTime(-1.5)
191  {}
192 
194  virtual ~StepperNewmarkImplicitAFormModifierXTest(){}
195 
197  virtual void modify(
199  const double time, const double dt,
201  {
202  switch(modType) {
204  {
205  testX_BEGIN_STEP = true;
206  testX = get_ele(*(x), 0);
207  break;
208  }
210  {
211  testX_BEFORE_SOLVE = true;
212  testDt = dt;
213  break;
214  }
216  {
217  testX_AFTER_SOLVE = true;
218  testTime = time;
219  testX = get_ele(*(x), 0);
220  break;
221  }
223  {
224  testX_END_STEP = true;
225  testX = get_ele(*(x), 0);
226  break;
227  }
228  default:
229  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
230  "Error - unknown action location.\n");
231  }
232  }
233 
234  bool testX_BEGIN_STEP;
235  bool testX_BEFORE_SOLVE;
236  bool testX_AFTER_SOLVE;
237  bool testX_END_STEP;
238  double testX;
239  double testXDot;
240  double testDt;
241  double testTime;
242 };
243 
244 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, AppAction_Modifier)
245 {
246  using Teuchos::RCP;
247  using Teuchos::getParametersFromXmlFile;
248  using Teuchos::sublist;
250 
251  double dt = 1.0;
252 
253  // Read params from .xml file
254  RCP<ParameterList> pList = getParametersFromXmlFile(
255  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
256  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
257 
258  // Setup the HarmonicOscillatorModel
259  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
262 
263  // Setup Stepper for field solve ----------------------------
265  Tempus::createStepperNewmarkImplicitAForm(model, Teuchos::null);
266 
267  auto modifier = rcp(new StepperNewmarkImplicitAFormModifierTest());
268  stepper->setAppAction(modifier);
269  stepper->initialize();
270 
271  // Setup TimeStepControl ------------------------------------
272  RCP<Tempus::TimeStepControl<double> > timeStepControl =
274  ParameterList tscPL = pl->sublist("Default Integrator")
275  .sublist("Time Step Control");
276  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
277  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
278  timeStepControl->setFinalTime(dt);
279  timeStepControl->setInitTimeStep(dt);
280  timeStepControl->initialize();
281 
282  // Setup initial condition SolutionState --------------------
283  using Teuchos::rcp_const_cast;
284  auto inArgsIC = model->getNominalValues();
286  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
288  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
289  RCP<Thyra::VectorBase<double> > icXDotDot =
290  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
292  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
293  icState->setTime (timeStepControl->getInitTime());
294  icState->setIndex (timeStepControl->getInitIndex());
295  icState->setTimeStep(0.0);
296  icState->setOrder (stepper->getOrder());
297  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
298 
299  // Setup SolutionHistory ------------------------------------
300  RCP<Tempus::SolutionHistory<double> > solutionHistory =
302  solutionHistory->setName("Forward States");
303  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
304  solutionHistory->setStorageLimit(2);
305  solutionHistory->addState(icState);
306 
307  // Setup Integrator -----------------------------------------
309  Tempus::createIntegratorBasic<double>();
310  integrator->setStepper(stepper);
311  integrator->setTimeStepControl(timeStepControl);
312  integrator->setSolutionHistory(solutionHistory);
313  integrator->initialize();
314 
315  // Integrate to timeMax
316  bool integratorStatus = integrator->advanceTime();
317  TEST_ASSERT(integratorStatus)
318 
319 
320  // Testing that each ACTION_LOCATION has been called.
321  TEST_COMPARE(modifier->testBEGIN_STEP, ==, true);
322  TEST_COMPARE(modifier->testBEFORE_SOLVE, ==, true);
323  TEST_COMPARE(modifier->testAFTER_SOLVE, ==, true);
324  TEST_COMPARE(modifier->testEND_STEP, ==, true);
325 
326  // Testing that values can be set through the Modifier.
327  auto x = integrator->getX();
328  auto Dt = integrator->getTime();
329  TEST_FLOATING_EQUALITY(modifier->testDt, Dt, 1.0e-14);
330  TEST_FLOATING_EQUALITY(modifier->testCurrentValue, get_ele(*(x), 0), 1.0e-14);
331  TEST_COMPARE(modifier->testName, ==, stepper->getStepperName());
332 }
333 
334 
335 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, AppAction_ModifierX)
336 {
337  using Teuchos::RCP;
338  using Teuchos::getParametersFromXmlFile;
339  using Teuchos::sublist;
341 
342  double dt = 1.0;
343 
344  // Read params from .xml file
345  RCP<ParameterList> pList = getParametersFromXmlFile(
346  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
347  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
348 
349  // Setup the HarmonicOscillatorModel
350  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
353 
354  // Setup Stepper for field solve ----------------------------
356  Tempus::createStepperNewmarkImplicitAForm(model, Teuchos::null);
357 
358  auto modifierX = rcp(new StepperNewmarkImplicitAFormModifierXTest());
359  stepper->setAppAction(modifierX);
360  stepper->initialize();
361 
362  // Setup TimeStepControl ------------------------------------
363  RCP<Tempus::TimeStepControl<double> > timeStepControl =
365  ParameterList tscPL = pl->sublist("Default Integrator")
366  .sublist("Time Step Control");
367  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
368  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
369  timeStepControl->setFinalTime(dt);
370  timeStepControl->setInitTimeStep(dt);
371  timeStepControl->initialize();
372 
373  // Setup initial condition SolutionState --------------------
374  using Teuchos::rcp_const_cast;
375  auto inArgsIC = model->getNominalValues();
377  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
379  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
380  RCP<Thyra::VectorBase<double> > icXDotDot =
381  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
383  Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
384  icState->setTime (timeStepControl->getInitTime());
385  icState->setIndex (timeStepControl->getInitIndex());
386  icState->setTimeStep(0.0);
387  icState->setOrder (stepper->getOrder());
388  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
389 
390  // Setup SolutionHistory ------------------------------------
391  RCP<Tempus::SolutionHistory<double> > solutionHistory =
393  solutionHistory->setName("Forward States");
394  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
395  solutionHistory->setStorageLimit(2);
396  solutionHistory->addState(icState);
397 
398  // Setup Integrator -----------------------------------------
400  Tempus::createIntegratorBasic<double>();
401  integrator->setStepper(stepper);
402  integrator->setTimeStepControl(timeStepControl);
403  integrator->setSolutionHistory(solutionHistory);
404  integrator->initialize();
405 
406  // Integrate to timeMax
407  bool integratorStatus = integrator->advanceTime();
408  TEST_ASSERT(integratorStatus)
409 
410  // Testing that each ACTION_LOCATION has been called.
411  TEST_COMPARE(modifierX->testX_BEGIN_STEP, ==, true);
412  TEST_COMPARE(modifierX->testX_BEFORE_SOLVE, ==, true);
413  TEST_COMPARE(modifierX->testX_AFTER_SOLVE, ==, true);
414  TEST_COMPARE(modifierX->testX_END_STEP, ==, true);
415 
416 
417  // Testing that values can be set through the Modifier.
418  auto Dt = integrator->getTime();
419  TEST_FLOATING_EQUALITY(modifierX->testDt, Dt, 1.0e-14);
420 
421  const auto x = integrator->getX();
422  TEST_FLOATING_EQUALITY(modifierX->testX, get_ele(*(x), 0), 1.0e-14);
423 }
424 
425 } // namespace Tempus_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)
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)