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