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