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