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