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