Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_HHTAlphaTest.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_StepperHHTAlpha.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::rcp;
41 using Teuchos::rcp_const_cast;
42 using Teuchos::ParameterList;
43 using Teuchos::sublist;
44 using Teuchos::getParametersFromXmlFile;
45 
46 using Tempus::IntegratorBasic;
47 using Tempus::SolutionHistory;
48 using Tempus::SolutionState;
49 
50 // Comment out any of the following tests to exclude from build/run.
51 #define TEST_BALL_PARABOLIC
52 #define TEST_CONSTRUCTING_FROM_DEFAULTS
53 #define TEST_SINCOS_SECONDORDER
54 #define TEST_SINCOS_FIRSTORDER
55 #define TEST_SINCOS_CD
56 
57 
58 #ifdef TEST_BALL_PARABOLIC
59 // ************************************************************
60 // ************************************************************
61 TEUCHOS_UNIT_TEST(HHTAlpha, BallParabolic)
62 {
63  //Tolerance to check if test passed
64  double tolerance = 1.0e-14;
65  // Read params from .xml file
66  RCP<ParameterList> pList =
67  getParametersFromXmlFile("Tempus_HHTAlpha_BallParabolic.xml");
68 
69  // Setup the HarmonicOscillatorModel
70  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
71  auto model = 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 
76  RCP<Tempus::IntegratorBasic<double> > integrator =
77  Tempus::integratorBasic<double>(pl, model);
78 
79  // Integrate to timeMax
80  bool integratorStatus = integrator->advanceTime();
81  TEST_ASSERT(integratorStatus)
82 
83  // Test if at 'Final Time'
84  double time = integrator->getTime();
85  double timeFinal =pl->sublist("Default Integrator")
86  .sublist("Time Step Control").get<double>("Final Time");
87  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
88 
89  // Time-integrated solution and the exact solution
90  RCP<Thyra::VectorBase<double> > x = integrator->getX();
91  RCP<const Thyra::VectorBase<double> > x_exact =
92  model->getExactSolution(time).get_x();
93 
94  // Plot sample solution and exact solution
95  std::ofstream ftmp("Tempus_HHTAlpha_BallParabolic.dat");
96  ftmp.precision(16);
97  RCP<const SolutionHistory<double> > solutionHistory =
98  integrator->getSolutionHistory();
99  bool passed = true;
100  double err = 0.0;
101  RCP<const Thyra::VectorBase<double> > x_exact_plot;
102  for (int i=0; i<solutionHistory->getNumStates(); i++) {
103  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
104  double time_i = solutionState->getTime();
105  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
106  x_exact_plot = model->getExactSolution(time_i).get_x();
107  ftmp << time_i << " "
108  << get_ele(*(x_plot), 0) << " "
109  << get_ele(*(x_exact_plot), 0) << std::endl;
110  if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
111  err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
112  }
113  ftmp.close();
114  std::cout << "Max error = " << err << "\n \n";
115  if (err > tolerance)
116  passed = false;
117 
118  TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
119  "\n Test failed! Max error = " << err << " > tolerance = " << tolerance << "\n!");
120 
121 }
122 #endif // TEST_BALL_PARABOLIC
123 
124 
125 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
126 // ************************************************************
127 // ************************************************************
128 TEUCHOS_UNIT_TEST(HHTAlpha, ConstructingFromDefaults)
129 {
130  double dt = 0.05;
131 
132  // Read params from .xml file
133  RCP<ParameterList> pList =
134  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_SecondOrder.xml");
135  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
136 
137  // Setup the HarmonicOscillatorModel
138  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
139  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
140 
141  // Setup Stepper for field solve ----------------------------
142  auto stepper = rcp(new Tempus::StepperHHTAlpha<double>());
143  stepper->setModel(model);
144  stepper->initialize();
145 
146  // Setup TimeStepControl ------------------------------------
147  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
148  ParameterList tscPL = pl->sublist("Default Integrator")
149  .sublist("Time Step Control");
150  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
151  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
152  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
153  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
154  timeStepControl->setInitTimeStep(dt);
155  timeStepControl->initialize();
156 
157  // Setup initial condition SolutionState --------------------
158  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
159  stepper->getModel()->getNominalValues();
160  auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
161  auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
162  auto icXDotDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
163  auto icState = rcp(new Tempus::SolutionState<double>(icX, icXDot, icXDotDot));
164  icState->setTime (timeStepControl->getInitTime());
165  icState->setIndex (timeStepControl->getInitIndex());
166  icState->setTimeStep(0.0);
167  icState->setOrder (stepper->getOrder());
168  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
169 
170  // Setup SolutionHistory ------------------------------------
171  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
172  solutionHistory->setName("Forward States");
173  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
174  solutionHistory->setStorageLimit(2);
175  solutionHistory->addState(icState);
176 
177  // Setup Integrator -----------------------------------------
178  RCP<Tempus::IntegratorBasic<double> > integrator =
179  Tempus::integratorBasic<double>();
180  integrator->setStepperWStepper(stepper);
181  integrator->setTimeStepControl(timeStepControl);
182  integrator->setSolutionHistory(solutionHistory);
183  //integrator->setObserver(...);
184  integrator->initialize();
185 
186 
187  // Integrate to timeMax
188  bool integratorStatus = integrator->advanceTime();
189  TEST_ASSERT(integratorStatus)
190 
191 
192  // Test if at 'Final Time'
193  double time = integrator->getTime();
194  double timeFinal =pl->sublist("Default Integrator")
195  .sublist("Time Step Control").get<double>("Final Time");
196  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
197 
198  // Time-integrated solution and the exact solution
199  RCP<Thyra::VectorBase<double> > x = integrator->getX();
200  RCP<const Thyra::VectorBase<double> > x_exact =
201  model->getExactSolution(time).get_x();
202 
203  // Calculate the error
204  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
205  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
206 
207  // Check the order and intercept
208  std::cout << " Stepper = " << stepper->description() << std::endl;
209  std::cout << " =========================" << std::endl;
210  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
211  std::cout << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
212  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
213  std::cout << " =========================" << std::endl;
214  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.144918, 1.0e-4 );
215 }
216 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
217 
218 
219 #ifdef TEST_SINCOS_SECONDORDER
220 // ************************************************************
221 TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_SecondOrder)
222 {
223  RCP<Tempus::IntegratorBasic<double> > integrator;
224  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
225  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
226  std::vector<double> StepSize;
227  std::vector<double> xErrorNorm;
228  std::vector<double> xDotErrorNorm;
229  const int nTimeStepSizes = 8;
230  double time = 0.0;
231 
232  // Read params from .xml file
233  RCP<ParameterList> pList =
234  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_SecondOrder.xml");
235 
236  // Setup the HarmonicOscillatorModel
237  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
238  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
239 
240  //Get k and m coefficients from model, needed for computing energy
241  double k = hom_pl->get<double>("x coeff k");
242  double m = hom_pl->get<double>("Mass coeff m");
243 
244  // Setup the Integrator and reset initial time step
245  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
246 
247  //Set initial time step = 2*dt specified in input file (for convergence study)
248  //
249  double dt =pl->sublist("Default Integrator")
250  .sublist("Time Step Control").get<double>("Initial Time Step");
251  dt *= 2.0;
252 
253  for (int n=0; n<nTimeStepSizes; n++) {
254 
255  //Perform time-step refinement
256  dt /= 2;
257  std::cout << "\n \n time step #" << n << " (out of "
258  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
259  pl->sublist("Default Integrator")
260  .sublist("Time Step Control").set("Initial Time Step", dt);
261  integrator = Tempus::integratorBasic<double>(pl, model);
262 
263  // Integrate to timeMax
264  bool integratorStatus = integrator->advanceTime();
265  TEST_ASSERT(integratorStatus)
266 
267  // Test if at 'Final Time'
268  time = integrator->getTime();
269  double timeFinal =pl->sublist("Default Integrator")
270  .sublist("Time Step Control").get<double>("Final Time");
271  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
272 
273  // Time-integrated solution and the exact solution
274  RCP<Thyra::VectorBase<double> > x = integrator->getX();
275  RCP<const Thyra::VectorBase<double> > x_exact =
276  model->getExactSolution(time).get_x();
277 
278  // Plot sample solution and exact solution at most-refined resolution
279  if (n == nTimeStepSizes-1) {
280  RCP<const SolutionHistory<double> > solutionHistory =
281  integrator->getSolutionHistory();
282  writeSolution("Tempus_HHTAlpha_SinCos_SecondOrder.dat", solutionHistory);
283 
284  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
285  for (int i=0; i<solutionHistory->getNumStates(); i++) {
286  double time_i = (*solutionHistory)[i]->getTime();
287  auto state = rcp(new Tempus::SolutionState<double>(
288  model->getExactSolution(time_i).get_x(),
289  model->getExactSolution(time_i).get_x_dot()));
290  state->setTime((*solutionHistory)[i]->getTime());
291  solnHistExact->addState(state);
292  }
293  writeSolution("Tempus_HHTAlpha_SinCos_SecondOrder-Ref.dat", solnHistExact);
294 
295  // Problem specific output
296  {
297  std::ofstream ftmp("Tempus_HHTAlpha_SinCos_SecondOrder-Energy.dat");
298  ftmp.precision(16);
299  RCP<const Thyra::VectorBase<double> > x_exact_plot;
300  for (int i=0; i<solutionHistory->getNumStates(); i++) {
301  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
302  double time_i = solutionState->getTime();
303  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
304  RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
305  x_exact_plot = model->getExactSolution(time_i).get_x();
306  //kinetic energy = 0.5*m*xdot*xdot
307  double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
308  ke *= 0.5*m;
309  //potential energy = 0.5*k*x*x
310  double pe = Thyra::dot(*x_plot, *x_plot);
311  pe *= 0.5*k;
312  double te = ke + pe;
313  //Output to file the following:
314  //[time, x computed, x exact, xdot computed, ke, pe, te]
315  ftmp << time_i << " "
316  << get_ele(*(x_plot), 0) << " "
317  << get_ele(*(x_exact_plot), 0) << " "
318  << get_ele(*(x_dot_plot), 0) << " "
319  << ke << " " << pe << " " << te << std::endl;
320  }
321  ftmp.close();
322  }
323  }
324 
325  // Store off the final solution and step size
326  StepSize.push_back(dt);
327  auto solution = Thyra::createMember(model->get_x_space());
328  Thyra::copy(*(integrator->getX()),solution.ptr());
329  solutions.push_back(solution);
330  auto solutionDot = Thyra::createMember(model->get_x_space());
331  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
332  solutionsDot.push_back(solutionDot);
333  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
334  StepSize.push_back(0.0);
335  auto solutionExact = Thyra::createMember(model->get_x_space());
336  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
337  solutions.push_back(solutionExact);
338  auto solutionDotExact = Thyra::createMember(model->get_x_space());
339  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
340  solutionDotExact.ptr());
341  solutionsDot.push_back(solutionDotExact);
342  }
343  }
344 
345  // Check the order and intercept
346  double xSlope = 0.0;
347  double xDotSlope = 0.0;
348  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
349  double order = stepper->getOrder();
350  writeOrderError("Tempus_HHTAlpha_SinCos_SecondOrder-Error.dat",
351  stepper, StepSize,
352  solutions, xErrorNorm, xSlope,
353  solutionsDot, xDotErrorNorm, xDotSlope);
354 
355  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
356  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00644755, 1.0e-4 );
357  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
358  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.104392, 1.0e-4 );
359 }
360 #endif // TEST_SINCOS_SECONDORDER
361 
362 
363 #ifdef TEST_SINCOS_FIRSTORDER
364 // ************************************************************
365 // ************************************************************
366 TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_FirstOrder)
367 {
368  RCP<Tempus::IntegratorBasic<double> > integrator;
369  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
370  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
371  std::vector<double> StepSize;
372  std::vector<double> xErrorNorm;
373  std::vector<double> xDotErrorNorm;
374  const int nTimeStepSizes = 8;
375  double time = 0.0;
376 
377  // Read params from .xml file
378  RCP<ParameterList> pList =
379  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_FirstOrder.xml");
380 
381  // Setup the HarmonicOscillatorModel
382  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
383  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
384 
385  //Get k and m coefficients from model, needed for computing energy
386  double k = hom_pl->get<double>("x coeff k");
387  double m = hom_pl->get<double>("Mass coeff m");
388 
389  // Setup the Integrator and reset initial time step
390  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
391 
392  //Set initial time step = 2*dt specified in input file (for convergence study)
393  //
394  double dt =pl->sublist("Default Integrator")
395  .sublist("Time Step Control").get<double>("Initial Time Step");
396  dt *= 2.0;
397 
398  for (int n=0; n<nTimeStepSizes; n++) {
399 
400  //Perform time-step refinement
401  dt /= 2;
402  std::cout << "\n \n time step #" << n << " (out of "
403  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
404  pl->sublist("Default Integrator")
405  .sublist("Time Step Control").set("Initial Time Step", dt);
406  integrator = Tempus::integratorBasic<double>(pl, model);
407 
408  // Integrate to timeMax
409  bool integratorStatus = integrator->advanceTime();
410  TEST_ASSERT(integratorStatus)
411 
412  // Test if at 'Final Time'
413  time = integrator->getTime();
414  double timeFinal =pl->sublist("Default Integrator")
415  .sublist("Time Step Control").get<double>("Final Time");
416  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
417 
418  // Time-integrated solution and the exact solution
419  RCP<Thyra::VectorBase<double> > x = integrator->getX();
420  RCP<const Thyra::VectorBase<double> > x_exact =
421  model->getExactSolution(time).get_x();
422 
423  // Plot sample solution and exact solution at most-refined resolution
424  if (n == nTimeStepSizes-1) {
425  RCP<const SolutionHistory<double> > solutionHistory =
426  integrator->getSolutionHistory();
427  writeSolution("Tempus_HHTAlpha_SinCos_FirstOrder.dat", solutionHistory);
428 
429  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
430  for (int i=0; i<solutionHistory->getNumStates(); i++) {
431  double time_i = (*solutionHistory)[i]->getTime();
432  auto state = rcp(new Tempus::SolutionState<double>(
433  model->getExactSolution(time_i).get_x(),
434  model->getExactSolution(time_i).get_x_dot()));
435  state->setTime((*solutionHistory)[i]->getTime());
436  solnHistExact->addState(state);
437  }
438  writeSolution("Tempus_HHTAlpha_SinCos_FirstOrder-Ref.dat", solnHistExact);
439 
440  // Problem specific output
441  {
442  std::ofstream ftmp("Tempus_HHTAlpha_SinCos_FirstOrder-Energy.dat");
443  ftmp.precision(16);
444  RCP<const Thyra::VectorBase<double> > x_exact_plot;
445  for (int i=0; i<solutionHistory->getNumStates(); i++) {
446  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
447  double time_i = solutionState->getTime();
448  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
449  RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
450  x_exact_plot = model->getExactSolution(time_i).get_x();
451  //kinetic energy = 0.5*m*xdot*xdot
452  double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
453  ke *= 0.5*m;
454  //potential energy = 0.5*k*x*x
455  double pe = Thyra::dot(*x_plot, *x_plot);
456  pe *= 0.5*k;
457  double te = ke + pe;
458  //Output to file the following:
459  //[time, x computed, x exact, xdot computed, ke, pe, te]
460  ftmp << time_i << " "
461  << get_ele(*(x_plot), 0) << " "
462  << get_ele(*(x_exact_plot), 0) << " "
463  << get_ele(*(x_dot_plot), 0) << " "
464  << ke << " " << pe << " " << te << std::endl;
465  }
466  ftmp.close();
467  }
468  }
469 
470  // Store off the final solution and step size
471  StepSize.push_back(dt);
472  auto solution = Thyra::createMember(model->get_x_space());
473  Thyra::copy(*(integrator->getX()),solution.ptr());
474  solutions.push_back(solution);
475  auto solutionDot = Thyra::createMember(model->get_x_space());
476  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
477  solutionsDot.push_back(solutionDot);
478  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
479  StepSize.push_back(0.0);
480  auto solutionExact = Thyra::createMember(model->get_x_space());
481  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
482  solutions.push_back(solutionExact);
483  auto solutionDotExact = Thyra::createMember(model->get_x_space());
484  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
485  solutionDotExact.ptr());
486  solutionsDot.push_back(solutionDotExact);
487  }
488  }
489 
490  // Check the order and intercept
491  double xSlope = 0.0;
492  double xDotSlope = 0.0;
493  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
494  double order = stepper->getOrder();
495  writeOrderError("Tempus_HHTAlpha_SinCos_FirstOrder-Error.dat",
496  stepper, StepSize,
497  solutions, xErrorNorm, xSlope,
498  solutionsDot, xDotErrorNorm, xDotSlope);
499 
500  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
501  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.048932, 1.0e-4 );
502  TEST_FLOATING_EQUALITY( xDotSlope, 1.18873, 0.01 );
503  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.393504, 1.0e-4 );
504 
505 }
506 #endif // TEST_SINCOS_FIRSTORDER
507 
508 
509 #ifdef TEST_SINCOS_CD
510 // ************************************************************
511 // ************************************************************
512 TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_CD)
513 {
514  RCP<Tempus::IntegratorBasic<double> > integrator;
515  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
516  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
517  std::vector<double> StepSize;
518  std::vector<double> xErrorNorm;
519  std::vector<double> xDotErrorNorm;
520  const int nTimeStepSizes = 8;
521  double time = 0.0;
522 
523  // Read params from .xml file
524  RCP<ParameterList> pList =
525  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_ExplicitCD.xml");
526 
527  // Setup the HarmonicOscillatorModel
528  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
529  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
530 
531  //Get k and m coefficients from model, needed for computing energy
532  double k = hom_pl->get<double>("x coeff k");
533  double m = hom_pl->get<double>("Mass coeff m");
534 
535  // Setup the Integrator and reset initial time step
536  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
537 
538  //Set initial time step = 2*dt specified in input file (for convergence study)
539  //
540  double dt =pl->sublist("Default Integrator")
541  .sublist("Time Step Control").get<double>("Initial Time Step");
542  dt *= 2.0;
543 
544  for (int n=0; n<nTimeStepSizes; n++) {
545 
546  //Perform time-step refinement
547  dt /= 2;
548  std::cout << "\n \n time step #" << n << " (out of "
549  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
550  pl->sublist("Default Integrator")
551  .sublist("Time Step Control").set("Initial Time Step", dt);
552  integrator = Tempus::integratorBasic<double>(pl, model);
553 
554  // Integrate to timeMax
555  bool integratorStatus = integrator->advanceTime();
556  TEST_ASSERT(integratorStatus)
557 
558  // Test if at 'Final Time'
559  time = integrator->getTime();
560  double timeFinal =pl->sublist("Default Integrator")
561  .sublist("Time Step Control").get<double>("Final Time");
562  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
563 
564  // Time-integrated solution and the exact solution
565  RCP<Thyra::VectorBase<double> > x = integrator->getX();
566  RCP<const Thyra::VectorBase<double> > x_exact =
567  model->getExactSolution(time).get_x();
568 
569  // Plot sample solution and exact solution at most-refined resolution
570  if (n == nTimeStepSizes-1) {
571  RCP<const SolutionHistory<double> > solutionHistory =
572  integrator->getSolutionHistory();
573  writeSolution("Tempus_HHTAlpha_SinCos_ExplicitCD.dat", solutionHistory);
574 
575  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
576  for (int i=0; i<solutionHistory->getNumStates(); i++) {
577  double time_i = (*solutionHistory)[i]->getTime();
578  auto state = rcp(new Tempus::SolutionState<double>(
579  model->getExactSolution(time_i).get_x(),
580  model->getExactSolution(time_i).get_x_dot()));
581  state->setTime((*solutionHistory)[i]->getTime());
582  solnHistExact->addState(state);
583  }
584  writeSolution("Tempus_HHTAlpha_SinCos_ExplicitCD-Ref.dat", solnHistExact);
585 
586  // Problem specific output
587  {
588  std::ofstream ftmp("Tempus_HHTAlpha_SinCos_ExplicitCD-Energy.dat");
589  ftmp.precision(16);
590  RCP<const Thyra::VectorBase<double> > x_exact_plot;
591  for (int i=0; i<solutionHistory->getNumStates(); i++) {
592  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
593  double time_i = solutionState->getTime();
594  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
595  RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
596  x_exact_plot = model->getExactSolution(time_i).get_x();
597  //kinetic energy = 0.5*m*xdot*xdot
598  double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
599  ke *= 0.5*m;
600  //potential energy = 0.5*k*x*x
601  double pe = Thyra::dot(*x_plot, *x_plot);
602  pe *= 0.5*k;
603  double te = ke + pe;
604  //Output to file the following:
605  //[time, x computed, x exact, xdot computed, ke, pe, te]
606  ftmp << time_i << " "
607  << get_ele(*(x_plot), 0) << " "
608  << get_ele(*(x_exact_plot), 0) << " "
609  << get_ele(*(x_dot_plot), 0) << " "
610  << ke << " " << pe << " " << te << std::endl;
611  }
612  ftmp.close();
613  }
614  }
615 
616  // Store off the final solution and step size
617  StepSize.push_back(dt);
618  auto solution = Thyra::createMember(model->get_x_space());
619  Thyra::copy(*(integrator->getX()),solution.ptr());
620  solutions.push_back(solution);
621  auto solutionDot = Thyra::createMember(model->get_x_space());
622  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
623  solutionsDot.push_back(solutionDot);
624  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
625  StepSize.push_back(0.0);
626  auto solutionExact = Thyra::createMember(model->get_x_space());
627  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
628  solutions.push_back(solutionExact);
629  auto solutionDotExact = Thyra::createMember(model->get_x_space());
630  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
631  solutionDotExact.ptr());
632  solutionsDot.push_back(solutionDotExact);
633  }
634  }
635 
636  // Check the order and intercept
637  double xSlope = 0.0;
638  double xDotSlope = 0.0;
639  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
640  double order = stepper->getOrder();
641  writeOrderError("Tempus_HHTAlpha_SinCos_ExplicitCD-Error.dat",
642  stepper, StepSize,
643  solutions, xErrorNorm, xSlope,
644  solutionsDot, xDotErrorNorm, xDotSlope);
645 
646  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
647  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00451069, 1.0e-4 );
648  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
649  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0551522, 1.0e-4 );
650 }
651 #endif // TEST_SINCOS_CD
652 
653 }
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