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