Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_NewmarkTest.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 
9 #include "Teuchos_UnitTestHarness.hpp"
10 #include "Teuchos_XMLParameterListHelpers.hpp"
11 #include "Teuchos_TimeMonitor.hpp"
12 
13 #include "Tempus_config.hpp"
14 #include "Tempus_IntegratorBasic.hpp"
15 #include "Tempus_StepperFactory.hpp"
16 
17 #include "../TestModels/HarmonicOscillatorModel.hpp"
18 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
19 
20 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
21 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
22 #include "Thyra_DetachedVectorView.hpp"
23 #include "Thyra_DetachedMultiVectorView.hpp"
24 
25 
26 #ifdef Tempus_ENABLE_MPI
27 #include "Epetra_MpiComm.h"
28 #else
29 #include "Epetra_SerialComm.h"
30 #endif
31 
32 #include <fstream>
33 #include <limits>
34 #include <sstream>
35 #include <vector>
36 
37 namespace Tempus_Test {
38 
39 using Teuchos::RCP;
40 using Teuchos::ParameterList;
41 using Teuchos::sublist;
42 using Teuchos::getParametersFromXmlFile;
43 
44 using Tempus::IntegratorBasic;
45 using Tempus::SolutionHistory;
46 using Tempus::SolutionState;
47 
48 //IKT, 3/22/17: comment out any of the following
49 //if you wish not to build/run all the test cases.
50 #define TEST_BALL_PARABOLIC
51 #define TEST_SINCOS_EXPLICIT
52 #define TEST_HARMONIC_OSCILLATOR_DAMPED_EXPLICIT
53 #define TEST_HARMONIC_OSCILLATOR_DAMPED_CTOR
54 #define TEST_HARMONIC_OSCILLATOR_DAMPED_SECOND_ORDER
55 #define TEST_HARMONIC_OSCILLATOR_DAMPED_FIRST_ORDER
56 
57 
58 #ifdef TEST_BALL_PARABOLIC
59 // ************************************************************
60 TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, BallParabolic)
61 {
62  //Tolerance to check if test passed
63  double tolerance = 1.0e-14;
64  // Read params from .xml file
65  RCP<ParameterList> pList =
66  getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_BallParabolic.xml");
67 
68  // Setup the HarmonicOscillatorModel
69  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
70  RCP<HarmonicOscillatorModel<double> > model =
71  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
72 
73  // Setup the Integrator and reset initial time step
74  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
75  RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
76  stepperPL->remove("Zero Initial Guess");
77 
78  RCP<Tempus::IntegratorBasic<double> > integrator =
79  Tempus::integratorBasic<double>(pl, model);
80 
81  // Integrate to timeMax
82  bool integratorStatus = integrator->advanceTime();
83  TEST_ASSERT(integratorStatus)
84 
85 // Test if at 'Final Time'
86  double time = integrator->getTime();
87  double timeFinal =pl->sublist("Default Integrator")
88  .sublist("Time Step Control").get<double>("Final Time");
89  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
90 
91  // Time-integrated solution and the exact solution
92  RCP<Thyra::VectorBase<double> > x = integrator->getX();
93  RCP<const Thyra::VectorBase<double> > x_exact =
94  model->getExactSolution(time).get_x();
95 
96  // Plot sample solution and exact solution
97  std::ofstream ftmp("Tempus_NewmarkExplicitAForm_BallParabolic.dat");
98  ftmp.precision(16);
99  RCP<const SolutionHistory<double> > solutionHistory =
100  integrator->getSolutionHistory();
101  bool passed = true;
102  double err = 0.0;
103  RCP<const Thyra::VectorBase<double> > x_exact_plot;
104  for (int i=0; i<solutionHistory->getNumStates(); i++) {
105  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
106  double time_i = solutionState->getTime();
107  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
108  x_exact_plot = model->getExactSolution(time_i).get_x();
109  ftmp << time_i << " "
110  << get_ele(*(x_plot), 0) << " "
111  << get_ele(*(x_exact_plot), 0) << std::endl;
112  if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
113  err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
114  }
115  ftmp.close();
116  std::cout << "Max error = " << err << "\n \n";
117  if (err > tolerance)
118  passed = false;
119 
120  TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
121  "\n Test failed! Max error = " << err << " > tolerance = " << tolerance << "\n!");
122 }
123 #endif
124 
125 
126 #ifdef TEST_SINCOS_EXPLICIT
127 // ************************************************************
128 TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, SinCos)
129 {
130  RCP<Tempus::IntegratorBasic<double> > integrator;
131  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
132  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
133  std::vector<double> StepSize;
134  std::vector<double> xErrorNorm;
135  std::vector<double> xDotErrorNorm;
136  const int nTimeStepSizes = 9;
137  double time = 0.0;
138 
139  // Read params from .xml file
140  RCP<ParameterList> pList =
141  getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_SinCos.xml");
142 
143  // Setup the HarmonicOscillatorModel
144  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
145  RCP<HarmonicOscillatorModel<double> > model =
146  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
147 
148 
149  // Setup the Integrator and reset initial time step
150  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
151  RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
152  stepperPL->remove("Zero Initial Guess");
153 
154  //Set initial time step = 2*dt specified in input file (for convergence study)
155  //
156  double dt =pl->sublist("Default Integrator")
157  .sublist("Time Step Control").get<double>("Initial Time Step");
158  dt *= 2.0;
159 
160  for (int n=0; n<nTimeStepSizes; n++) {
161 
162  //Perform time-step refinement
163  dt /= 2;
164  std::cout << "\n \n time step #" << n << " (out of "
165  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
166  pl->sublist("Default Integrator")
167  .sublist("Time Step Control").set("Initial Time Step", dt);
168  integrator = Tempus::integratorBasic<double>(pl, model);
169 
170  // Integrate to timeMax
171  bool integratorStatus = integrator->advanceTime();
172  TEST_ASSERT(integratorStatus)
173 
174  // Test if at 'Final Time'
175  time = integrator->getTime();
176  double timeFinal =pl->sublist("Default Integrator")
177  .sublist("Time Step Control").get<double>("Final Time");
178  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
179 
180  // Plot sample solution and exact solution
181  if (n == 0) {
182  RCP<const SolutionHistory<double> > solutionHistory =
183  integrator->getSolutionHistory();
184  writeSolution("Tempus_NewmarkExplicitAForm_SinCos.dat", solutionHistory);
185 
186  RCP<Tempus::SolutionHistory<double> > solnHistExact =
187  Teuchos::rcp(new Tempus::SolutionHistory<double>());
188  for (int i=0; i<solutionHistory->getNumStates(); i++) {
189  double time_i = (*solutionHistory)[i]->getTime();
190  RCP<Tempus::SolutionState<double> > state =
191  Teuchos::rcp(new Tempus::SolutionState<double>(
192  model->getExactSolution(time_i).get_x(),
193  model->getExactSolution(time_i).get_x_dot()));
194  state->setTime((*solutionHistory)[i]->getTime());
195  solnHistExact->addState(state);
196  }
197  writeSolution("Tempus_NewmarkExplicitAForm_SinCos-Ref.dat",solnHistExact);
198  }
199 
200  // Store off the final solution and step size
201  StepSize.push_back(dt);
202  auto solution = Thyra::createMember(model->get_x_space());
203  Thyra::copy(*(integrator->getX()),solution.ptr());
204  solutions.push_back(solution);
205  auto solutionDot = Thyra::createMember(model->get_x_space());
206  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
207  solutionsDot.push_back(solutionDot);
208  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
209  StepSize.push_back(0.0);
210  auto solutionExact = Thyra::createMember(model->get_x_space());
211  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
212  solutions.push_back(solutionExact);
213  auto solutionDotExact = Thyra::createMember(model->get_x_space());
214  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
215  solutionDotExact.ptr());
216  solutionsDot.push_back(solutionDotExact);
217  }
218  }
219 
220  // Check the order and intercept
221  double xSlope = 0.0;
222  double xDotSlope = 0.0;
223  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
224  double order = stepper->getOrder();
225  writeOrderError("Tempus_NewmarkExplicitAForm_SinCos-Error.dat",
226  stepper, StepSize,
227  solutions, xErrorNorm, xSlope,
228  solutionsDot, xDotErrorNorm, xDotSlope);
229 
230  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
231  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
232  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
233  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.233045, 1.0e-4 );
234 
235  Teuchos::TimeMonitor::summarize();
236 }
237 #endif
238 
239 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED_EXPLICIT
240 // ************************************************************
241 TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, HarmonicOscillatorDamped)
242 {
243  RCP<Tempus::IntegratorBasic<double> > integrator;
244  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
245  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
246  std::vector<double> StepSize;
247  std::vector<double> xErrorNorm;
248  std::vector<double> xDotErrorNorm;
249  const int nTimeStepSizes = 9;
250  double time = 0.0;
251 
252  // Read params from .xml file
253  RCP<ParameterList> pList =
254  getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml");
255 
256  // Setup the HarmonicOscillatorModel
257  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
258  RCP<HarmonicOscillatorModel<double> > model =
259  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
260 
261 
262  // Setup the Integrator and reset initial time step
263  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
264  RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
265  stepperPL->remove("Zero Initial Guess");
266 
267  //Set initial time step = 2*dt specified in input file (for convergence study)
268  //
269  double dt =pl->sublist("Default Integrator")
270  .sublist("Time Step Control").get<double>("Initial Time Step");
271  dt *= 2.0;
272 
273  for (int n=0; n<nTimeStepSizes; n++) {
274 
275  //Perform time-step refinement
276  dt /= 2;
277  std::cout << "\n \n time step #" << n << " (out of "
278  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
279  pl->sublist("Default Integrator")
280  .sublist("Time Step Control").set("Initial Time Step", dt);
281  integrator = Tempus::integratorBasic<double>(pl, model);
282 
283  // Integrate to timeMax
284  bool integratorStatus = integrator->advanceTime();
285  TEST_ASSERT(integratorStatus)
286 
287  // Test if at 'Final Time'
288  time = integrator->getTime();
289  double timeFinal =pl->sublist("Default Integrator")
290  .sublist("Time Step Control").get<double>("Final Time");
291  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
292 
293  // Plot sample solution and exact solution
294  if (n == 0) {
295  RCP<const SolutionHistory<double> > solutionHistory =
296  integrator->getSolutionHistory();
297  writeSolution("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.dat", solutionHistory);
298 
299  RCP<Tempus::SolutionHistory<double> > solnHistExact =
300  Teuchos::rcp(new Tempus::SolutionHistory<double>());
301  for (int i=0; i<solutionHistory->getNumStates(); i++) {
302  double time_i = (*solutionHistory)[i]->getTime();
303  RCP<Tempus::SolutionState<double> > state =
304  Teuchos::rcp(new Tempus::SolutionState<double>(
305  model->getExactSolution(time_i).get_x(),
306  model->getExactSolution(time_i).get_x_dot()));
307  state->setTime((*solutionHistory)[i]->getTime());
308  solnHistExact->addState(state);
309  }
310  writeSolution("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Ref.dat", solnHistExact);
311  }
312 
313  // Store off the final solution and step size
314  StepSize.push_back(dt);
315  auto solution = Thyra::createMember(model->get_x_space());
316  Thyra::copy(*(integrator->getX()),solution.ptr());
317  solutions.push_back(solution);
318  auto solutionDot = Thyra::createMember(model->get_x_space());
319  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
320  solutionsDot.push_back(solutionDot);
321  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
322  StepSize.push_back(0.0);
323  auto solutionExact = Thyra::createMember(model->get_x_space());
324  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
325  solutions.push_back(solutionExact);
326  auto solutionDotExact = Thyra::createMember(model->get_x_space());
327  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
328  solutionDotExact.ptr());
329  solutionsDot.push_back(solutionDotExact);
330  }
331  }
332 
333  // Check the order and intercept
334  double xSlope = 0.0;
335  double xDotSlope = 0.0;
336  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
337  //double order = stepper->getOrder();
338  writeOrderError("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Error.dat",
339  stepper, StepSize,
340  solutions, xErrorNorm, xSlope,
341  solutionsDot, xDotErrorNorm, xDotSlope);
342 
343  TEST_FLOATING_EQUALITY( xSlope, 1.060930, 0.01 );
344  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.508229, 1.0e-4 );
345  TEST_FLOATING_EQUALITY( xDotSlope, 1.019300, 0.01 );
346  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.172900, 1.0e-4 );
347 
348  Teuchos::TimeMonitor::summarize();
349 }
350 #endif
351 
352 
353 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED_CTOR
354 // ************************************************************
355 // ************************************************************
356 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, ConstructingFromDefaults)
357 {
358  double dt = 1.0;
359 
360  // Read params from .xml file
361  RCP<ParameterList> pList = getParametersFromXmlFile(
362  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
363  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
364 
365  // Setup the HarmonicOscillatorModel
366  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
367  RCP<HarmonicOscillatorModel<double> > model =
368  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
369 
370  // Setup Stepper for field solve ----------------------------
371  auto sf = Teuchos::rcp(new Tempus::StepperFactory<double>());
372  RCP<Tempus::StepperNewmarkImplicitAForm<double> > stepper =
373  sf->createStepperNewmarkImplicitAForm(model, Teuchos::null);
374 
375  // Setup TimeStepControl ------------------------------------
376  RCP<Tempus::TimeStepControl<double> > timeStepControl =
377  Teuchos::rcp(new Tempus::TimeStepControl<double>());
378  ParameterList tscPL = pl->sublist("Default Integrator")
379  .sublist("Time Step Control");
380  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
381  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
382  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
383  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
384  timeStepControl->setInitTimeStep(dt);
385  timeStepControl->initialize();
386 
387  // Setup initial condition SolutionState --------------------
388  using Teuchos::rcp_const_cast;
389  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
390  stepper->getModel()->getNominalValues();
391  RCP<Thyra::VectorBase<double> > icX =
392  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
393  RCP<Thyra::VectorBase<double> > icXDot =
394  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
395  RCP<Thyra::VectorBase<double> > icXDotDot =
396  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
397  RCP<Tempus::SolutionState<double> > icState =
398  Teuchos::rcp(new Tempus::SolutionState<double>(icX, icXDot, icXDotDot));
399  icState->setTime (timeStepControl->getInitTime());
400  icState->setIndex (timeStepControl->getInitIndex());
401  icState->setTimeStep(0.0);
402  icState->setOrder (stepper->getOrder());
403  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
404 
405  // Setup SolutionHistory ------------------------------------
406  RCP<Tempus::SolutionHistory<double> > solutionHistory =
407  Teuchos::rcp(new Tempus::SolutionHistory<double>());
408  solutionHistory->setName("Forward States");
409  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
410  solutionHistory->setStorageLimit(2);
411  solutionHistory->addState(icState);
412 
413  // Setup Integrator -----------------------------------------
414  RCP<Tempus::IntegratorBasic<double> > integrator =
415  Tempus::integratorBasic<double>();
416  integrator->setStepperWStepper(stepper);
417  integrator->setTimeStepControl(timeStepControl);
418  integrator->setSolutionHistory(solutionHistory);
419  //integrator->setObserver(...);
420  integrator->initialize();
421 
422 
423  // Integrate to timeMax
424  bool integratorStatus = integrator->advanceTime();
425  TEST_ASSERT(integratorStatus)
426 
427 
428  // Test if at 'Final Time'
429  double time = integrator->getTime();
430  double timeFinal =pl->sublist("Default Integrator")
431  .sublist("Time Step Control").get<double>("Final Time");
432  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
433 
434  // Time-integrated solution and the exact solution
435  RCP<Thyra::VectorBase<double> > x = integrator->getX();
436  RCP<const Thyra::VectorBase<double> > x_exact =
437  model->getExactSolution(time).get_x();
438 
439  // Calculate the error
440  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
441  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
442 
443  // Check the order and intercept
444  std::cout << " Stepper = " << stepper->description() << std::endl;
445  std::cout << " =========================" << std::endl;
446  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
447  std::cout << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
448  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
449  std::cout << " =========================" << std::endl;
450  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4 );
451 }
452 #endif
453 
454 
455 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED_SECOND_ORDER
456 // ************************************************************
457 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, HarmonicOscillatorDamped_SecondOrder)
458 {
459  RCP<Tempus::IntegratorBasic<double> > integrator;
460  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
461  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
462  std::vector<double> StepSize;
463  std::vector<double> xErrorNorm;
464  std::vector<double> xDotErrorNorm;
465  const int nTimeStepSizes = 10;
466  double time = 0.0;
467 
468  // Read params from .xml file
469  RCP<ParameterList> pList =
470  getParametersFromXmlFile("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
471 
472  // Setup the HarmonicOscillatorModel
473  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
474  RCP<HarmonicOscillatorModel<double> > model =
475  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
476 
477 
478  // Setup the Integrator and reset initial time step
479  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
480 
481  //Set initial time step = 2*dt specified in input file (for convergence study)
482  //
483  double dt =pl->sublist("Default Integrator")
484  .sublist("Time Step Control").get<double>("Initial Time Step");
485  dt *= 2.0;
486 
487  for (int n=0; n<nTimeStepSizes; n++) {
488 
489  //Perform time-step refinement
490  dt /= 2;
491  std::cout << "\n \n time step #" << n << " (out of "
492  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
493  pl->sublist("Default Integrator")
494  .sublist("Time Step Control").set("Initial Time Step", dt);
495  integrator = Tempus::integratorBasic<double>(pl, model);
496 
497  // Integrate to timeMax
498  bool integratorStatus = integrator->advanceTime();
499  TEST_ASSERT(integratorStatus)
500 
501  // Test if at 'Final Time'
502  time = integrator->getTime();
503  double timeFinal =pl->sublist("Default Integrator")
504  .sublist("Time Step Control").get<double>("Final Time");
505  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
506 
507  // Plot sample solution and exact solution
508  if (n == 0) {
509  RCP<const SolutionHistory<double> > solutionHistory =
510  integrator->getSolutionHistory();
511  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.dat", solutionHistory);
512 
513  RCP<Tempus::SolutionHistory<double> > solnHistExact =
514  Teuchos::rcp(new Tempus::SolutionHistory<double>());
515  for (int i=0; i<solutionHistory->getNumStates(); i++) {
516  double time_i = (*solutionHistory)[i]->getTime();
517  RCP<Tempus::SolutionState<double> > state =
518  Teuchos::rcp(new Tempus::SolutionState<double>(
519  model->getExactSolution(time_i).get_x(),
520  model->getExactSolution(time_i).get_x_dot()));
521  state->setTime((*solutionHistory)[i]->getTime());
522  solnHistExact->addState(state);
523  }
524  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder-Ref.dat", solnHistExact);
525  }
526 
527  // Store off the final solution and step size
528  StepSize.push_back(dt);
529  auto solution = Thyra::createMember(model->get_x_space());
530  Thyra::copy(*(integrator->getX()),solution.ptr());
531  solutions.push_back(solution);
532  auto solutionDot = Thyra::createMember(model->get_x_space());
533  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
534  solutionsDot.push_back(solutionDot);
535  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
536  StepSize.push_back(0.0);
537  auto solutionExact = Thyra::createMember(model->get_x_space());
538  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
539  solutions.push_back(solutionExact);
540  auto solutionDotExact = Thyra::createMember(model->get_x_space());
541  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
542  solutionDotExact.ptr());
543  solutionsDot.push_back(solutionDotExact);
544  }
545  }
546 
547  // Check the order and intercept
548  double xSlope = 0.0;
549  double xDotSlope = 0.0;
550  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
551  double order = stepper->getOrder();
552  writeOrderError("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder_SinCos-Error.dat",
553  stepper, StepSize,
554  solutions, xErrorNorm, xSlope,
555  solutionsDot, xDotErrorNorm, xDotSlope);
556 
557  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
558  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0484483, 1.0e-4 );
559  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
560  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0484483, 1.0e-4 );
561 
562  Teuchos::TimeMonitor::summarize();
563 }
564 #endif
565 
566 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED_FIRST_ORDER
567 // ************************************************************
568 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, HarmonicOscillatorDamped_FirstOrder)
569 {
570  RCP<Tempus::IntegratorBasic<double> > integrator;
571  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
572  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
573  std::vector<double> StepSize;
574  std::vector<double> xErrorNorm;
575  std::vector<double> xDotErrorNorm;
576  const int nTimeStepSizes = 10;
577  double time = 0.0;
578 
579  // Read params from .xml file
580  RCP<ParameterList> pList =
581  getParametersFromXmlFile("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.xml");
582 
583  // Setup the HarmonicOscillatorModel
584  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
585  RCP<HarmonicOscillatorModel<double> > model =
586  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
587 
588 
589  // Setup the Integrator and reset initial time step
590  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
591 
592  //Set initial time step = 2*dt specified in input file (for convergence study)
593  //
594  double dt =pl->sublist("Default Integrator")
595  .sublist("Time Step Control").get<double>("Initial Time Step");
596  dt *= 2.0;
597 
598  for (int n=0; n<nTimeStepSizes; n++) {
599 
600  //Perform time-step refinement
601  dt /= 2;
602  std::cout << "\n \n time step #" << n << " (out of "
603  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
604  pl->sublist("Default Integrator")
605  .sublist("Time Step Control").set("Initial Time Step", dt);
606  integrator = Tempus::integratorBasic<double>(pl, model);
607 
608  // Integrate to timeMax
609  bool integratorStatus = integrator->advanceTime();
610  TEST_ASSERT(integratorStatus)
611 
612  // Test if at 'Final Time'
613  time = integrator->getTime();
614  double timeFinal =pl->sublist("Default Integrator")
615  .sublist("Time Step Control").get<double>("Final Time");
616  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
617 
618  // Plot sample solution and exact solution
619  if (n == 0) {
620  RCP<const SolutionHistory<double> > solutionHistory =
621  integrator->getSolutionHistory();
622  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.dat", solutionHistory);
623 
624  RCP<Tempus::SolutionHistory<double> > solnHistExact =
625  Teuchos::rcp(new Tempus::SolutionHistory<double>());
626  for (int i=0; i<solutionHistory->getNumStates(); i++) {
627  double time_i = (*solutionHistory)[i]->getTime();
628  RCP<Tempus::SolutionState<double> > state =
629  Teuchos::rcp(new Tempus::SolutionState<double>(
630  model->getExactSolution(time_i).get_x(),
631  model->getExactSolution(time_i).get_x_dot()));
632  state->setTime((*solutionHistory)[i]->getTime());
633  solnHistExact->addState(state);
634  }
635  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Ref.dat", solnHistExact);
636  }
637 
638  // Store off the final solution and step size
639  StepSize.push_back(dt);
640  auto solution = Thyra::createMember(model->get_x_space());
641  Thyra::copy(*(integrator->getX()),solution.ptr());
642  solutions.push_back(solution);
643  auto solutionDot = Thyra::createMember(model->get_x_space());
644  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
645  solutionsDot.push_back(solutionDot);
646  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
647  StepSize.push_back(0.0);
648  auto solutionExact = Thyra::createMember(model->get_x_space());
649  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
650  solutions.push_back(solutionExact);
651  auto solutionDotExact = Thyra::createMember(model->get_x_space());
652  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
653  solutionDotExact.ptr());
654  solutionsDot.push_back(solutionDotExact);
655  }
656  }
657 
658  // Check the order and intercept
659  double xSlope = 0.0;
660  double xDotSlope = 0.0;
661  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
662  double order = stepper->getOrder();
663  writeOrderError("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Error.dat",
664  stepper, StepSize,
665  solutions, xErrorNorm, xSlope,
666  solutionsDot, xDotErrorNorm, xDotSlope);
667 
668  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
669  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0224726, 1.0e-4 );
670  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
671  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0122223, 1.0e-4 );
672 
673  Teuchos::TimeMonitor::summarize();
674 }
675 #endif
676 }
Tempus_StepperFactory.hpp
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:306
Tempus_Test::HarmonicOscillatorModel
Consider the ODE: where is a constant, is a constant damping parameter, is a constant forcing par...
Definition: HarmonicOscillatorModel_decl.hpp:52
Tempus_Test::writeOrderError
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar > > stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope)
Definition: Tempus_ConvergenceTestUtils.hpp:184
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Definition: Tempus_Integrator.hpp:26
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:512
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:30
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56